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FOREWORD 


The  requirement  for  designing  geophysical  surveys  In  such  a way  that 
the  true  continuous  field  may  be  derived  to  a predetermined  accuracy  from 
digital  samples  Is  basic  to  any  data  collection  operation.  The  time 
required  In  the  collection  and  analysis  of  the  reconnaissance  data,  which 
Is  necessary  for  the  utilization  of  these  survey  design  procedures,  may  be 
Justified  not  only  from  a purely  scientific  standpoint  but,  with  the 
escalating  cost  of  shipboard  operations,  from  an  economic  standpoint  as 
well.  Application  of  the  theory  and  algorithms  presented  In  this  report 
Is  expected  to  result  In  a significant  Improvement  In  survey  platform 
utilization. 


J.  E.  AWffiS 
Captain,  U.S.  Navy 
Commander 
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ABSTRACT 

A theory  for  designing  parallel  track-type  geophysical  surveys,  as  • 
well  as  the  necessary  numerical  algorithms  for  Implementing  this 
theory,  has  been  developed  which  is  easily  applied  to  many  diff^irent 
sampling  orobisnij.  Within  this  context,  survey  design  consists  of 
defining  the  appropriate  track  spacing  track  direction,  and  down-track 
sampling  rate  which  will  produce  a set  of  discrete  digital  measure- 
ments describing  the  environment  to  a predetermined  accuracy. 

These  procedures  are  based  primarily  upon  the  use  of  one  and 
two-dimensional  Fourier  transforms  applied  to  appropriate  numerical 
models  of  the  sampling  process  in  order  tc  estimate  the  variance  or 
mean  square  erroi  as  well  as  the  spectral  content  of  the  sampling 
error.  Since  thc.se  error  estimates  are  computed  in  the  spatial 
frequency  domain,  application  of  the  convolution  theorem  is  shown  to 
produce  a particularly  efficient  process  for  propagating  the  error  es- 
timates through  a variety  of  linear  operations  performed  upon  the 
survey  data. 

Several  practical  applications  are  presented  to  illustrate  the 
adaptability  .of  the  theory.  These  applications  include  the  ncf'r  real- 
time design  of  hydrographic  surveys  utilizing  a small-scale  computer, 
the  design  of  gravity  surveys  from  which  estimates  of  vertical  deflec- 
tion and  geoid  undulation  may  be  derived  to  a specified  accuracy,  and 
the  design  of  oceanic  sound  speed  surveys  which  illustrates  the 
eppliCbtion  of  the  theory  to  three-dimensional  fields. 
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INTRODUCTION 


In  any  area  of  science  which  relies  upon  the  utilization  of 
discrete  measurements  of  a spatially  or  temporally  variable  physical 
process  as  a primary  source  of  information,  it  is  highly  desirable,  if 
not  critical,  that  estimates  be  made  of  the  error  which  is  introduced 
by  discrete  sampling.  In  many  cases,  a requirement  exists  for  esti- 
mating not  only  this  primary  sampling  error  associated  with  the 
original  measurements,  but  a secondary  sampling  error  which  is  the 
result  of  propagating  the  primary  error  through  a series  of  linear 
operations  to  produce  quantities  computed  from  the  sampled  data. 

Examples  of  this  type  of  linear  operation  are  upward  and  downward  ana- 
lytic continuation  of  potential  fields,  one  and  two-dimensional 
differentiation,  and  the  calculation  of  vertical  deflection  components 
and  geoidal  undulations  from  measurements  of  the  earth's  gravity 
field. 

Although  the  appropriate  estimates  of  sampling  error  may  be  comput- 
ed after  the  measurements  are  completed,  by  far  the  most  efficient  use 
of  the  theory  and  algorithms  described  in  this  report  is  made  during 
the  design  phase  of  the  data  collection  operation.  The  application  of 
these  procedures  to  geophysical  processes  falls  under  the  general 
heading  of  geophysical  survey  design.  Large  scale  geophysical  surveys 
are  presently  conducted  primarily  by  aircraft,  ships,  and  satellites 
travelling  along  survey  tracks  laid  out  in  some  prescribed  manner. 

Under  these  conditions,  geophysical  survey  design  consists  of  defining 
the  appropriate  track  spacing,  track  direction,  and  down-track  sampling 
rate  which  will  produce  a set  of  discrete  measurements  describing  the 

environment  to  a predetermined  accuracy. 
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Some  preliminary  theory  for  estimating  sampling  error  has  been 
developed  by  applying  statistical  methods  based  on  one-dimensional 
covariance  functions  to  physical  geodesy.  These  methods  were  utilized 
by  Rapp  (1964)  to  estimate  numerically  the  error  of  prediction  or 
Interpolation  of  gravity  anomalies.  Heiskanen  and  Moritz  (1967)  out- 
lined the  procedure  for  propagating  this  covariance  estimate  of  sam- 
pling error  through  the  calculation  of  spherical  harmonics  of  the 
gravity  field  and  the  computation  of  geoidal  undulations.  The  assump- 
tion of  isotropy  inherent  in  the  application  of  one-dimensional 
covariance  functions  to  two-dimensional  fieltls,  the  practical  problems 
of  computing  accurate  local  covariance  estimates  utilizing  a finite 
amount  of  data  from  a process  which  is  basically  non-stationary  with 
respect  to  spatial  coordinates,  and  the  numerical  complexity  of 
error  propagation  through  linear  operations  severely  limits  the 
utility  of  these  existing  procedures  for  survey  design. 

In  order  to  overcome  these  limitations,  a theory  of  geophysical 
survey  design,  as  well  as  the  necessary  numerical  algorithm  for  imple- 
menting this  theory,  has  been  developed  which  is  quite  easily  applied 
to  many  different  sampling  problems.  This  theory  is  particularly 
efficient  for  propagating  the  desired  sampling  error  estimates  through 
a variety  of  linear  operations.  These  procedures  are  based  primarily 
upon  the  use  of  one  and  two-dimensional  Fourier  transforms  applied  to 
.’umerical  models  of  the  sampling  process  in  order  to  estimate  the 
variance  or  mean  square  error  as  well  as  the  spectral  content  of  the 
sampling  error.  Since  these  error  estimates  are  computed  in  the 
spatial  frequency  domain,  straightforward  application  of  the  conv'olu- 
tion  theorem  results  in  an  efficient  error  propagation  process. 

Several  practical  applications  are  presented  to  demonstrate  the 
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use  of  the  theory  and  the  relative  ease  with  which  it  may  be  adapted 
to  many  different  types  of  survey  operations.  Although  in  some  cases, 
sampling  error  may  be  considered  to  have  a minor  impact  upon  the  quality 
of  the  results  obtained  from  an  analysis  of  the  data,  this  does  not 
negate  the  requirement  for  designing  an  efficient  survey  plan,  especial- 
ly in  light  of  the  present  cost  of  large  scale  survey  operations. 
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THEORY  OF  SURVEY  DESIGN 


In  order  to  develop  a comprehensive  theory  which  could  form  the 
basis  for  a practical  procedure  to  design  geophysical  surveys,  the 
following  general  criteria  were  used: 

1.  The  theory  should  be  easily  adaptable  to  a wide  variety 
of  survey  operations,  with  the  emphasis  being  placed  on 
the  design  of  large  scale  track-type  geophysical 
surveys. 

2.  Application  of  the  theory  should  allow  estimates  to  be 
made  of  the  total  variance  or  mean  square  error  and  the 
two-dimensional  spectral  content  of  the  sampling  error 
as  a function  of  the  survey  pattern. 

3.  The  error  estimates  should  be  in  a form  that  would 
facilitate  their  propagation  through  a variety  of 
linear  operations. 

4.  The  computational  burden  associated  with  the  numerical 
implementation  of  the  theory  should  be  minimized  to  the 
extent  that  survey  design  could  be  performed  in  near 
real-time  with  small  scale  computers. 

Utilizing  these  general  criteria  as  guidelines,  a procedure  for 
’ 'ing  track-type  surveys  has  been  developed  which  is  based  upon  the 
use  of  mcthcmatical  models  of  the  survey  process,  combined  with  classi- 
cal sampling  theory,  to  estimate  the  two-dimensional  spectrum  of  the 
error  cS*-oci with  the  discrete  sampling  of  a continuous  two-  or 
three-dimensional  function  which  is  essentially  non-stationary.  The 
successful  practical  implementation  of  this  theory  involves  the  use  of 

reconnaissance  data  to  estimate  the  boundaries  of  statistically 
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homogeneous  provinces,  the  design  of  appropriate  models  of  the  survey 
pattern,  and  the  reliable  estimation  of  both  one-  and  two-dimensional 
Fourier  transforms  of  geophysical  data.  In  this  section,  the  theory, 
and  the  numerical  formulation  developed  for  implementation  will  be  pre- 
sented in  detail.  Fortran  IV  computer  programs  are  available  from  the 
author  for  each  of  the  algorithms  required  for  practical  utilization 
of  this  theory. 

Practical  Problems  in  Spectral  Estimation 

Throughout  this  discussion,  the  following  definitions  will  be  used. 

The  direct  Fourier  transform  is  defined  by  G(to)  = g(x)e  dx, 

1.  itiix 

and  the  inverse  Fourier  transform  is  defined  by  g(x)  = ® dw. 

The  problems  associated  with  estimating  the  one-dimensional  auto- 
correlation function  and  power  spectrum  of  a continuous  function  from  a 
sample  of  finite  length  are  presented  in  detail  in  the  classical  work 
by  Blackman  and  Tukey  (1958).  More  recently,  Lacoss  (1971)  applied  two 
new  nonlinear  spectral  analysis  techniques  to  the  problem  of  resolving 
narrow  spectral  peaks  in  seismic  data.  These  techniques  (Maximum 
Likelihood  Method  developed  by  Capon  (1969) , and  Maximum  Entropy 
Method  utilized  by  Burg  (1967)),  appear  to  yield  results  superior  to 
those  obtained  by  straightfoni/ard  window  modification  applied  to  the 
estimated  autocorrelation  function. 

In  contrast,  the  problem  of  estimating  the  spectral  content  of  a 
finite  length  of  data  from  geophysical  processes  such  as  the  earth's 
gravity  or  magnetic  fields  is  generated  not  by  prominent  spectral 
peaks  but  bj’  the  large  amount  of  energy  contained  in  wavelengths  con- 
siderably longer  than  the  length  of  the  available  data  record. 
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Processes  which  possess  this  type  of  spectral  character  have  been  term 
cd  "red  noise"  processes  by  Shapiro  and  Ward  (1960) . The  generaliza- 
tion of  a pure  "red  noise"  spectrum  as  the  power  spectrum  of  a first- 
order  linear  Markov  process  has  been  developed  by  Gilman  et  al  (1963) . 
Alldrcdge,  et  al  (1963)  computed  a one-dimensional  spectrum  of  the 
earth's  magnetic  field  which  shows  the  large  amount  of  energy  in  the 
field  at  the  long  wavelengths.  In  a statistical  sense,  this  long  wave 
length  energy  may  be  considered  as  the  non-stationary  component  of  our 
finite  length  sample.  Since  the  basis  of  survey  design  lies  in  esti- 
mating the  Fourier  transform  of  the  error  generated  by  sampling  a 
locally  stationary  two-  or  three-dimensional  function»  two  distinct 
problems  must  be  addressed.  These  problems  are  outlined  as  follows: 

1.  A significant  amount  of  the  energy  contained  in  the  spectra  of 
the  earth's  gravity  and  magnetic  fields  is  located  in  wave- 
lengths much  longer  than  the  extent  of  any  particular  survey 
operation.  This  long  - wavelength  energy  causes  g(x,y)  to 
appear  to  be  non-stationary  with  regard  to  spatial  coordinates, 
that  is,  the  mean  and  autocorrelation  function  of  g(x,y)  vary 
over  the  survey  area. 

2.  Only  in  rare  instances  are  closely  spaced  gridded  values  of 
g(x,y)  available  for  numerically  computing  a local  two- 
dimensional  spectrum  which  may  be  assumed  to  be  free  of 
saraplirr  error. 

In  practice,  the  second  problem  is  generally  annoying  in  that  it 

requires  tiiu  collection  of  additiona_  survey  data,  but  the  first 

problem  is  often  fatal  unless  adequate  precautions  are  taken.  The 

first  problem  generates  what  is  often  referred  to  as  spectral  leakage 

(Tukey,  1967).  In  one-dimension,  the  collection  of  a finite  length  of 
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data  (-T/2  T/2)  from  a continuous  function  r(x)  consists  essential- 

ly of  multiplying  g(x)  by  a rectangular  data  window  function  w(x) 
defined  by 


w(x) 


1,  -T/2  <x£  T/2 
0,  elsewhere 


It  may  be  shovm  by  straightforward  integration  that  the  Fourier  trans- 
form of  w(x)  is  given  by 


VJ(w) 


Sin  T cj/2 
T(o/2  ■’ 


where  T is  the  total  length  of  the  data  window. 

In  the  frequency  domain,  the  effect  of  multiplying  g(x)  by  w(x)  to 

obtain  a finite  length  of  data  is  equivalent  to  convolving  G(co)  with 

W(o)).  Thus,  the  Fourier  transform  which  is  actually  computed  is  not 

G(o))  but  G (oj) , which  is  given  by 
c 

G^(o))  = j G(u)  W((o-u)du. 


Historically,  many  different  data  window  functions  have  been  proposed 

and  used  in  an  attempt  to  make  G (w)  a better  approximation  of  G(w) 

c 

(see  e.g.  Blackman  and  Tukey,  1958).  Prior  to  the  current  trend  of 
computing  spectra  through  utilization  of  the  FFT  (Fast  Fourier  Trans- 
form) developed  by  Cooley  and  Tukey  (1965),  these  data  windows  were 
applied  to  the  sample  autocorrelation  function.  With  the  development 
of  the  FFT  as  a standard  analysis  tool,  the  data  windows  must  be 
applied  directly  to  the  sampled  function.  The  design  of  these  data 
window  functions  involves  a compromise  between  the  "sharpness"  of  the 
frequency  cutoff  of  W(oj)  and  the  magnitude  of  the  side  lobes.  For 
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example,  the  window  function  consisting  of  a so  called  cosine  taper 
defined  as 


w^(x)  = 


l/2(l+cos  ^),  - T/2  <x<T/2 
0,  elsewhere 


is  often  used  to  reduce  the  magnitude  of  the  side  lobes  at  the  expense 
of  broadening  the  main  lobe.  The  Fourier  transform  of  the  cosine  taper 
window  is  given  by 


„ = I Sin_Tu)/j  ^ T 

”t^“^  2 Tcj/2  ^ 4 


Sin  (tt-Tw/2)  . Sin  (7r+Tw/2) 
z: — r:: + n: — rz 


7T-Tw/2 


tt+Tw/2 


W(f)  ”t^^> 

Figure  1 is  a plot  of  • ' and  — ^ — in  which  the  substitution  w=27rf 

has  been  made  to  more  clearly  indicate  the  relationships  between  the 

transforms  of  the  window  functions  and  their  length  (T) . Note  that  the 

^t^^^  W(f) 

main  lobe  of  — — is  much  wider  than  that  of  — but  the  magnitude 

of  the  first  side  lobe  has  been  decreased  from  approximately  0.21  to 

approximately  0.01.  In  practice,  the  use  of  w^(x)  tends  to  reduce  the 

oscillatory  character  of  G (w) , but  a significant  amount  of  spectral 

c 

leakage  may  still  be  present.  As  a consequence  of  the  negative  side 
lobe  on  the  cosine  taper  window,  negative  power  spectral  estimates 
often  occur  when  this  window  is  applied  to  the  autocorrelation  function 
computed  from  processes  exhibiting  "red  noise"  character.  Other 
windows,  such  s'  *'he  Bartlett  win ’ow  (Blackman  and  Tukey,  1958),  defined 
as 

2|xl 


WgCx) 


1 - 


. Ixl<lT/2j, 


0,  elsewhere 


may  be  used  to  eliminate  the  negative  power  estimates  at  the  expense  of 

Increasing  the  magnitude  of  the  side  lobes. 
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Figure  1.  Fourier  transform  of  cosine  toper  (w^(t)A) 
and  rectangular  (w(f)A)  data  windows. 
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For  a real,  asymmetric  g(x),  the  Fourier  transform  is  a complex 

quantity  of  the  form  G(w)  = E(u)  + iD(w)  with  an  even  real  part  (E(w)) 

and  an  odd  imaginary  part  (D(w)).  The  amplitude  spectrum  of  g(x)  is 

defined  as  |g(oj)|  and  the  power  spectrum  or  energy  spectrum  is  defined 
2 

as  |g(w)|  (Bracewell,  1965).  The  computed  amplitude  spectrum  of  g(x) 
multiplied  by  a data  window  is  then  given  by 


|G  (u))1  = ([E(w)*W(u)]2  + [D0a))*W(w)]2)^''^ 

' c 

where  * denotes  the  convolution  operation.  In  the  case  of  gravity  and 
magnetic  fields,  this  leakage  of  the  long  wavelength  energy  through  the 
side  lobes  of  the  window  function  may  easily  become  serious  enough  to 
distort  the  computed  spectrum  completely.  In  light  of  the  ever  increas- 
ing utilization  of  the  FFT,  the  importance  of  this  problem  cannot  be 
overemphasized . 

As  a quantitative  illustration  of  this  leakage  effect,  a compari- 
son is  made  between  the  theoretical  and  numerical  amplitude  spectrum  of 
the  gravity  anomaly  generated  by  a deep-seated  two-dimensional  fault 
model  shoxjn  in  Figure  2.  The  gravity  anomaly  along  a profile  perpendic- 
ular to  the  strike  of  this  body  for  a dip  angle  (A)  of  90®  is  given  by 
'Davj,  1971) 


g(x)  ^ + 2yp|(D+T)  tan“'(-^)  - D tan"' ^)  + ^ In  (P^^)  j\ 


where  D,  x.  and  p are  defined  as  in  Figure  2,  and  y is  the  universal 

gravitatioiioj.  constant.  Since  g(x)  is  not  absolutely  integrable 

( lg(x)ldx  is  not  finite),  the  Fourier  transform  of  g(x)  theoretical- 
J — OO 

ly  does  not  exist.  Nevertheless,  a representation  of  the  transform 

which  will  be  sufficient  for  our  needs  is  given  by 
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Figure  2.  The  two-dimensional  gravity  fault  model. 
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|g(co)|  = 2ifYPS  (irT6  (w) ) ^+ 


1 , -WD  -U)(D+T). 
“Tie  -e  J 

0)2  ■' 


(Davis,  1971),  where 


6(co)  is  the  Dirac  delta  function.  For  w/o,  this  reduces  to 


!<=(-)  I 


-0)D  -o)(D+T) 

e -e 


The  curve  denoted  by  the  symbol  (A)  in 
Figure  3 is  a plot  of  the  theoretical  amplitude  spectrum  of  the  gravity 
anomaly  generated  by  a fault  model  for  a depth  of  1 km,  a thickness  of 
10  km,  and  a density  contrast  of  1 g/cc.  Digital  data  for  computing  a 
numerical  FFT  amplitude  spectrum  were  generated  by  evaluating  g(x)  at  a 
data  interval  (DI)  of  AX  = 0.5  km  for  the  range  -64  km  <x<  64  km.  Note 
that  all  FFT  spectra  are  computed  and  plotted  at  equal  increments  of 
normalized  frequency  with  units  of  cycles  per  data  interval  (*"/DI)  . 

The  resulting  numerical  spectrum  computed  from  this  data  set  using  the 
rectangular  window  function,  w(x),  is  shown  by  the  curve  denoted  by  the 
symbol  (o)  in  Figure  3.  The  effect,  on  the  computed  amplitude  spectrum, 
of  multiplying  g(x)  by  w(x)  with  the  resulting  convolution  in  the 
frequency  domain  is  readily  apparent.  The  distortion  of  the  numerical 
spectrum  at  the  high  frequency  end  is  generated  by  leakage  as  the  side 
lobes  of  W(u))  operate  on  those  spectral  components  of  G(o))  with  wave- 
lengths longer  than  128  km.  This  distortion  is  typical  of  that  obtain- 
ed in  practice  from  geophysical  processes  such  as  gravity  or  magnetic 
^ip-ds  which  possess  a "red  noise"  type  of  spectrum. 

In  practice,  an  alternative  procedure  called  prev;hitening  (Tukey, 
1967)  has  been  found  to  yield  more  consistently  accurate  spectral  esti- 
mates than  ui.ose  generally  obtairable  by  applying  window  modifications 
to  this  type  of  "red  noise"  process.  In  this  application,  prewhitening 
is  considered  to  be  a numerical  process  which  operates  on  the  original 


data  set  in  such  a way  that  the  resulting  amplitude  spectrum  is  nearly 
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NORMALIZED  FREQUENCY  (-^/Ol) 

Figure  3,  Numerical  example  of  spectral  leakage. 


flat  (within  one  or  two  orders  of  magnitude  for  gravity  or  magnetic 
data)  over  the  frequency  range  of  Interest  and  nearly  zero  outside 
this  range.  The  manner  in  which  the  prewhitening  operation  reduces 
the  effect  of  applying  a rectangular  data  window  may  be  shown  as 
follows ; 

Define  a function  g(x)  which  possesses  a flat  spectrum,  i.e., 

g(x)  = 6(x),  and  Giui)  = 1. 

With  0)  = 2irf»  G(f)  = 1,  and  the  Fourier  transform  of  the  rectangular 
data  window  is  given  by 


W(f) 


„ SinirfT 

T 

irfT 


The  computed  Fourier  transform  (G  (f))  Is  the  result  of  convolving 

c 

G(f)  with  W(f),  thus, 

0 (f)  . I r Sln£nl_ 
c ^ 'o  n 

and  (Gradshteyn  and  Ryzhlk,  1965) , 


G (f)  = G(f)  = 1. 
c 

As  is  shown  by  this  example,  the  application  of  the  rectangular  data 
window  to  a function  possessing  a flat  spectrum  has  no  adverse  effect 
upon  the  computed  spectrum. 

Assuming  that  N equally  spaced  digital  data  values  are  available 
for  analysis,  this  prewhitening  process  is  most  easily  accomplished  by 
applying  an  appropriately  designed  digital  filter  to  these  N data 
values,  oxuce,  as  the  example  shown  in  Figure  3 indicates,  a 
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realistic  spectrum  cannot  generally  be  computed  prior  to  prewhitening, 
the  appropriate  amplitude  response  of  this  filter  can  be  obtained  only 
through  experience  with  a particular  type  of  data.  As  a practical 
matter,  a technique  which  has  been  found  to  be  particularly  effective 
in  this  regard  is  to  compute  the  spectrum  of  a sample  set  of  data  using 
the  cosine  taper  window  w^(x)  and  then  develop  the  appropriate  response 
for  the  prew’hitening  filter  from  this  first  estimate.  Another  equally 
effective,  though  slightly  more  time  consuming  method,  consists  of 
applying  a set  of  contiguous  numerical  band-pass  filters  to  the  data 
set,  and  integrating  the  output  of  each  filter  to  produce  an  estimate 
of  the  energy  contained  within  each  frequency  band.  In  addition,  a 
priori  knowledge  of  the  expected  form  of  the  spectrum  for  the  type  of 
data  being  analyzed  is  sometimes  available.  In  the  case  of  two- 
dimensional  gravity  or  magnetic  fields,  the  fact  that  the  amplitude 
spectrum  is  generally  dominated  by  a term  of  the  form  e where  D is 
the  depth  to  the  disturbing  body,  is  helpful. 

A filter  possessing  an  amplitude  response  as  shown  in  Figure  3, 
which  has  been  set  to  effectively  remove  normalized  frequencies  less 
than  1/N  cycles  per  data  interval  has  been  found  to  work  reasonably 
well  for  both  gravity  and  magnetic  data.  Details  concerning  the  design 
of  this  type  of  filter  are  presented  in  Appendix  A.  Figure  3 illus- 
trates the  improvement  v;hich  was  obtained  by  prewhitening  the  gravity 
fault  model  data  in  the  space  domain,  computing  the  amplitude  spectrum 
of  these  filtered  data,  then  utilizing  the  convolution  theorem  to 
correct  this  spectrum  for  the  effect  of  the  filter  over  the  frequency 
range  l/N  to  0.5  cycles  per  data  intei*val  by  dividing  the  computed  FFT 

by  the  frequency  response  of  the  filter.  The  final  corrected  FFT 
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spectrum  (denoted  by  the  s>Tnbol  (x)  is  nearly  identical  to  the  theoret- 
ical curve  over  the  entire  frequency  range. 

The  FFT  spectrum  obtained  through  application  of  the  cosine  taper 
window  (w^(x))  is  also  presented  for  comparison.  Note  that  there  is 
essentially  no  residual  oscillation  in  the  low  frequency  components  of 
the  prewhitened  and  corrected  FFT  while  the  side  lobes  of  the  cosine 
taper  window  operating  on  the  6(w)  component  of  G(w)  produced  a signif- 
icant ripple  at  these  frequencies,  .^s  this  example  illustrates,  the 
design  of  the  prewhitening  filter  is  not  particularly  critical  (ie.  the 
prewhitened  spectrum  does  not  need  to  be  exactly  white),  but  the 
employment  of  some  form  of  prewhitening  or  window  modification  is 
imperative  if  accurate  spectral  estimates  are  to  be  expected  from  a 
finite  length  of  data  obtained  from  a "red  noise"  process. 

One-Dimensional  Sampling 

As  developed  here,  the  underlying  theory  utilized  in  geophysical 
survey  design  consists  of  combining  the  Shannon  sampling  theorem 
(Shannon,  1949)  with  various  models  of  one-  and  two-dimensional  sam- 
pling functions  and  Parseval's  relationship  (Hamming,  1962)  to  estimate 
the  variance  or  mean  square  error  and  the  spectral  content  of  the 
sampling  error  as  a function  of  survey  track  direction  and  track 
spacing. 

The  proof  of  the  one-dimensional  Shannon  sampling  theorem  is  out- 
lined ■’*"  ®apoulis  (1962).  Since  this  is  the  basis  for  practical  survey 
design,  the  detailed  proof  is  presented  in  Appendix  B.  This  theorem 
states  that,  given  equally  spaced  digital  values  f(mT)  of  a continuous 

function  f(t)  which  is  band-limited  to  the  normalized  frequency  range 
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• iT 


-Ti<to<ir,  and  such  that 


fTf 

|f(w)1 

^-TT 


dw  exists,  f(t)  may  be  represented  by 


^ .r/  Sin  n(t-raT)  , . 

f(t)  . ? tfaT)  ■■  “11  1- 

m=~oo  ' 


In  this  formulation,  T is  the  spacing  between  the  digital  samples  of 
f(t).  Thus,  this  discrete  numerical  convolution  of  the  digital  samples 
of  f(t)  with  a Sine  function  ( — amounts  to  a perfect  interpola- 
tion procedure. 

Within  the  context  of  survey  design,  the  ultimate  aim  is  to  esti- 
mate the  variance  and  spectral  content  of  the  error  generated  by 
sampling  a continuous  function,  v;hich  is  band-limited  for  practical 
purposes,  at  a rate  less  than  that  required  by  the  Shannon  theorem. 

In  order  to  generate  this  estimate  of  the  so-called  sampling  error, 
consider  that  the  equally  spaced  digital  samples  (f(mT))  of  the 
continuous  function  (f(t))  are  obtained  by  multiplying  f(t)  by  an 
infinite  sequence  of  unit  impulses  spaced  a distance  T apart. 

Bracewell  (1965)  called  this  sequence  the  shav  symbol  and  defines  it  as 

oo 

lll(t)  = Z o(t-mT).  Hsu  (1970)  shows,  in  a straightforv/ard  manner, 
that  the  Fourier  transform  of  lll(t)  = ^ Z 6(u)  - . 

m=-<» 

As  a consequence  of  the  convolution  theorem  for  Fourier  transforms 

(see  e.g.  Papoulis,  1962),  the  multiplication  of  f(t)  by  lll(t)  to 

produce  the  digital  samples  f(mT)  in  the  time  domain  results  in 

convolving  the  Fourier  transform  of  lll(t)  with  the  Fourier  transform 

of  f(t).  These  relationships  are  shown  in  Figure  4.  From  this  figure, 

it  is  clear  that,  as  T increases,  the  spacing  between  the  impulses  of 

lll(w)  becomes  less.  If  T is  allowed  to  increase  to  the  point  where 
2tt 

— <2A,  then  the  convolution  operation  replicates  F(w)  in  such  a way 
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Figure  4.  One-dimensional  sampling. 
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TT 

that  F(io)  overlaps  itself  at  the  folding  frequency  (±^)  producing  what 
is  called  sampling  error. 

The  Fourier  transform  of  this  sampling  error  is  composed  of  two 
parts.  The  first  part  may  be  considered  to  be  the  error  of  commission 
or  the  aliasing  error.  The  Fourier  transform  of  this  error  component 
is  defined  as  the  difference  between  F(w)  and  111  (u)  * F(w)  up  to  the 
folding  frequency.  The  amplitude  spectrum  of  this  error  component  is 
given  by 

|F^(w)|  = I 111(a))  * F(u)  - F(o))l  (1) 

The  second  part  of  the  sampling  error  may  be  considered  to  be  the 
error  of  omission  since  spectral  components  of  F(o))  for  ^ are 

not  recoverable  from  the  sampled  data.  The  amplitude  spectrum  of  this 
error  component  is  simply 

1f^(0))1  = |f(oj)|  for  y a. 

These  two  error  components  combine  to  form  the  sampling  error  spectrum 
defined  as 

|Fs(o))  I = |f^(o))  I + |f^(u))|  (2) 

This  error  is  shouTi  in  Figure  5 for  the  real  part  of  F(u). 

For  survey  design  operations,  the  interest  lies  not  only  in  the 
spectrum  of  the  sampling  error  but  also  in  the  mean  square  error  or 
RMS  of  the  error  as  a function  of  sample  spacing.  This  mean  square 
error  estimate  is  readily  available  througli  the  application  of  Parseval* 
formula  (Papoulis,  1962).  This  formula  states  that  energy  is  conserved 
when  transformations  are  made  between  the  time  or  space  domain  and  the 
frequency  domain  and  is  given  by 
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Figure  5.  The  two  components  of  one-dimensionol  sampling  error. 
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|G(f)  1 df  for  w = 2TTf  and 

^-A 


g(t)  bandlimited  to  ±A,  In  practice,  G(f)  is  computed  numerically  uti- 
lizing the  Fast  Fourier  transform  which  produces  discrete  estimates  of 
G(f)  at  a normalized  frequency  spacing  of  cycles  per  data  interval  for 
N digital  values  (g^)  of  g(t).  A numerical  estimate  of  the  mean  square 
sampling  error  using  Parseval's  formula  is  then  given  by 


1 

N 


N-1 

Z 

m=o 


m 


N/2-1 
+ 2 I 


K=1 


(3) 


where  G_(— ) are  the  numerical  FFT  estimates  of  the  sampling  error  spec- 
b N 

trum,  defined  by  equation  (2),  /or  J = 0,....N/2. 

With  regard  to  the  computation  of  the  amplitude  or  pov;er  spectrum 
of  the  aliasing  error  through  the  use  of  equation  (1),  an  approximation 
to  |f  (w)|  is  presented  here  for  use  in  real  time  survey  design  with 
small  scale  computers.  This  technique  is  particularly  appropriate  when 
only  the  sampling  error  variance  is  desired,  rather  than  the  sampling 
error  spectrum.  This  approximation  is  called  a power  fold,  and  is 
developed  in  the  following  manner. 

For  a real  asymmetric  f (x) , the  Fourier  transform  is  a complex 
quantity  of  the  form  F(u)  = E(oj)  + i D(to)  with  an  even,  real  part 
(E((J)),  and  an  odd,  imaginary  part  (D(w)).  A complex  quantity  posses- 
sing this  type  of  symmetry  is  termed  Hermitian.  An  explicit  formula- 
tion of  the  power  spectrum  of  the  aliasing  error  component  for  a 
sample  spacing  (T)  is,  from  equation  (1), 
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E(u)  L 6[(o)-u)- -^)du 
00  m=-«o 


|lll(io)*F(a))-F(a))l^  = ||2.  I 

I ' — B< 

#00 

+ 1 D(u)  ? 6[((0-u)-  ^]du  - [E((0)+iD(0))] 

J —00  ui=— 00 


Changing  the  order  of  integration  and  summation,  and  utilizing  the 

r 

property  of  6(a))  given  by  E (o)) 5 (o)) du)  = E(o)  , yields 

j —CO 

2 


|lll(a))*F(a))-F(a)) 


2iT  ~ 27rm. 

— Z E(0)  - — ) 

m=-°° 


2tt  ^ 

^ Z D(co- 

21TI11v 
T ' 

■ - D(0)) 

k 

• - E(o)) 


(A) 


In  the  power  fold  approximation,  the  computer  storage  requirement  and 

computation  time  is  essentially  cut  in  half  by  approximatljng  the  actual 

power  spectrum  of  the  aliasing  error  component  given  by  equation  (A) 

2 

by  replicating  |f(o))|  rather  than  F(cj)  . 

2 2 2 

Thus,  |lll(a))*F(a))-F(a))  I =^|lll(a))*|F(a))  ] - 1f(o))I  | for 

Explicitly,  we  have 


|111(U))*  |f(0))|^  - 1f(0))|^|* 


2^  ? r.2/-  2Trm. 

— Z EUdi-  — ) 

nj=s-co 


- E^(a)) 


I.  J 2tt  * ^2  / 2-rnn. 

Mt- 

m=-oo 


- D^(a))J- 


(5) 


Comparing  equations  (A)  and  (5),  it  is  readily  apparent  that  the 
accuracy  of  this  power  fold  approximation  is  a function  of  the  ampli- 
tude of  the  cross  terms  in  the  product  of  the  summations  in  equation 
(A).  Empirical  tests  on  many  different  types  of  geophysical  data 
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indicate  that  the  error  of  this  approximation  is  within  a small  percent- 
age of  the  total  mean  square  sampling  error  except  at  folding  frequen- 
cies equivalent  to  the  first  two  or  three  harmonics  of  the  sample  data 
record.  Two  examples  of  these  tests  are  shown  in  Figures  6 and  7. 

Figure  6 indicates  the  sampling  error  estimates  obtained  from  applying 
equations  (A)  and  (5)  to  a 32  point  numerical  test  function  shown  in 
Figure  8.  In  this  test,  the  power  fold  is  seen  to  yield  relatively 
accurate  estimates  of  the  sampling  error  except  for  sample  spacings 
associated  with  the  first  three  harmonics  of  the  32  point  test  function. 
Figure  7 shows  a similar  result  obtained  by  applying  equations  (A) 
and  (5)  to  a 128  point  digital  record  of  temporal  variations  of  oceanic 
sound  speed.  As  in  the  previous  example,  estimates  of  the  sampling 
error  for  sample  spacings  associated  with  folding  frequencies  equivalent 
to  the  first  three  harmonics  of  the  128  point  data  record  are  relative- 
ly Inaccurate. 

In  order  to  Illustrate  these  concepts  in  a numerical  example,  a 
digital  test  function  containing  predominantly  high  frequency  energy 
was  generated  utilizing  equations  (A-1)  and  (A-2)  from  Appendix  A to 
produce  a set  of  high-pass  filter  weights.  The  filter  control  para- 
meters were  set  to  V = 0.08,  H = 0.2,  N = 15.  The  numerical  weight 

c 

function  generated  by  this  process  is  sho\'m  in  Figure  8.  Note  that  the 

central  values  of  this  function  should  be  scaled  by  the  factor  100. 

2 

The  power  spectrum  of  this  function  [F(w) | , computed  via  the  FFT,  is 

shown  in  Figure  9 as  well  as  the  estimate  of  the  aliased  spectrum 

produced  for  a AX  sample  spacing  of  two  data  intervals.  This  estimate 

2 

was  produced  by  a numerical  solution  of  llll(u)*F(w) | for  the  real 

2 

and  imaginary  fold  and  lll(u)* |F(u) | for  the  power  fold  approxima- 
tion. In  this  example,  with  AX  = two  data  intervals,  the  folding 
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SAMIM.E  SFACiNG  (DATA  INTERVAIS) 

Figure  6.  Num«ncal  on«>dimensiofiol  test  - Sampling  •rror  vs.  sample  spocing. 
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SAMPLE  SPACING  (DATA  INTERVALS) 

Figure  7.  Example  of  sompling  error  for  temporal  voriotions  of  oceanic  sound  speed  at  a 
depth  of  75  meters — one  data  interval =25  min. 
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Figure  8.  One-dimensional  digital  test  function. 
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Figure  9.  Numerical  one-dimensional  test— Aliased  spectrum  for  aX=2DI. 
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frequency  is  0.25  cycles  per  data  interval.  Figure  10  is  a plot  of  the 
estimated  power  spectrum  of  the  sampling  error  obtained  by  a numerical 
solution  of  equations  (4)  and  (5).  At  this  sample  spacing,  the  sam- 
pling error  spectrum  produced  by  the  power  fold  approximation  was 
identical  to  that  produced  by  folding  the  real  and  imaginary  parts  of 
F(o)). 

As  a test  of  this  numerical  process,  the  original  test  function 
was  sampled  at  every  second  data  point  and  a cubic  spline  interpolation 
(Davis  and  Kontis,  1970)  was  applied  to  obtain  values  at  the  missing 
points.  As  shown  in  Figure  10,  the  power  spectrum  of  the  difference 
between  the  interpolated  values  and  the  true  values  of  the  test  func- 
tion is  in  excellent  agreement  with  that  predicted  by  the  theory.  The 
spectral  content  of  the  error  for  frequencies  less  than  0.25  cycles  per 
data  interval  is  the  error  of  commission  v:hile  the  error  for  frequencies 
greater  than  0.25  cycles  per  data  interval  is  the  error  of  omission. 

Utilizing  equation  (3) , the  estimated  mean  square  sampling  error 
for  both  the  real  and  imaginary  fold  and  the  power  fold  approximation 
is  0.021201,  The  actual  sampling  error  computed  from  the  interpolated 
data  values  was  0.019330.  In  practice,  the  sampling  error  spectrum  and 
mean  square  error  may  be  estimated  for  any  folding  frequency  which  is  a 
harmonic  of  the  original  data  record.  Figure  6 is  a plot  of  the  sam- 
pling error  estimates  for  all  possible  folding  frequencies  of  this 
test  function.  It  is  interesting  to  note  that,  with^the  total  sampling 
error  detxned  as  consisting  of  both  the  error  of  commission  and  the 
error  of  omission,  it  is  perfectly  feasible  to  generate  a mean  square 
sampling  error  that  is  greater  than  the  mean  square  value  of  the 

original  function.  In  this  particular  example,  the  mean  square  value 
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Figure  10.  Numerical  one  •dimensional  test— Spectrum  of  sampling  error  for  aX=2DI. 
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of  the  original  test  function  is  0.010625,  Referring  to  Figure  6,  it 
Is  clear  that  the  sampling  error  may  be  expected  to  exceed  this  value 
for  any  sample  spacing  greater  than  I.IA  data  intervals. 

Two-Dimensional  Track  Sampling 

At  the  present  time,  most  large  scale  geophysical  survey  opera- 
tions are  conducted  by  collecting  data  along  nominally  parallel  survey 
tracks  by  the  use  of  aircraft,  ships,  helicopters,  and  satellites. 

The  basic  concept  of  geophysical  survey  design  Involves  determining 
the  optimum  track  spacing,  track  direction,  and  do\<m-track  sample  rate 
which  will  produce  digital  data  to  describe  these  continuous  fields  to 
a predetermined  accuracy. 

With  the  current  widespread  use  of  digital  recording  systems  and 
the  utilization  of  the  one-dimensional  sampling  theory  for  selecting  the 
appropriate  down-track  sample  rate,  we  can  reasonably  assume  that  the 
down- track  sample  rate  is  sufficient  to  define  exactly  a two-dimensional 
function  f(x,y)  along  each  survey  track.  Under  this  assumption,  the 
appropriate  mathematical  model  for  defining  two-dimensional  track-type 
surveys  is  the  raster  sampling  model  utilized  by  Bracewell  (1965 ) as  a 
model  for  the  formation  of  television  images.  This  model  is  essential- 
ly a set  of  parallel  delta  function  ridges  generated  by  considering  the 
one-dimensional  shaw  symbol  as  a two-dimensional  function,  ie: 

CO 

s (x,y)  = g(x)  E 6(y-mT)  where  g(x)=l  for 

ni=-oo 

survey  tracks  parallel  to  the  x axis  and  spaced  a distance  T apart  in 
the  y direction. 

Since  the  Fourier  transform  of  a two-dimensional  function  which  is 
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a product  of  two  one-dimensional  functions  is  the  product  of  the  one- 


dimensional transforms,  the  Fourier  transform  of  S^(x,y)  is  given  by 


S^(u,v)  = 2tt6(u)  E 6(v-^) 


where  2ti6(u)  is  the  Fourier  transform  of  the  constant  g(x)=l.  This 
model  is  sho\im  in  Figure  11.  The  survey  data  D(x,y)  collected  along  a 
set  of  equally  spaced  survey  tracks  oriented  parallel  to  the  x axis  is 
then  given  by  D(x,y)  = f(x,y)  S^(x,y).  Since  multiplication  in  the 
two-dimensional  space  domain  results  in  a convolution  in  the  two- 
dimensional  spatial  frequency  domain,  the  result  of  the  survey  opera- 
tion is  to  replicate  the  true  two-dimensional  Fourier  transform  of 
f(x,y)  in  the  cross-track  direction  by  convolving  F(u,v)  with  S^(u,v)). 
As  in  the  one-dimensional  case,  the  two-dimensional  sampling  error  is 
defined  as  consisting  of  an  aliasing  error  or  error  of  commission  for 
frequencies  less  than  the  folding  frequency  and  an  error  of  omission  for 
frequencies  greater  than  the  folding  frequency.  The  real  part  of  the 
spectrum  of  this  sampling  error  (Fg(u,v))is  shown  in  Figure  11. 

An  explicit  formulation  for  the  two-dimensional  pov;er  spectrum  is 
given  by 


2 


2 


+ 


2 


(6) 
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where  R(u,v)  and  I(u,v)  arc  the  real  and  iraaginary  parts  of  F(u,v).  As 

in  the  one-dimensional  case,  the  error  of  omission  for  tracks  parallel 

2 2 

to  the  X axis  is  simply  [f^Cu.v)]  = |f(u,v)|  for  v greater  than  the 

folding  frequency.  It  should  be  noted  that,  in  contrast  to  the  one- 
dimensional  transform,  the  real  and  imaginary  parts  of  the  two-dimen- 
sional transform  of  an  arbitrary  f(x,y),  which  does  not  exhibit 
circular  symmetry,  are  composed  of  both  even  and  odd  components.  In 
practice,  the  numerical  evaluation  of  equation  (6)  is  faciliated  by 
separating  both  R(u,v)  and  I(u,v)  into  their  respective  even  and  odd 
components  by  use  of  the  identity  f(t)  = •|•(f  (t)+f  (-t)  ] for  the  even  com- 
ponent and  f(c)  = ^[f (t)-f (-t) ] for  the  odd  component  (Wylie,  1960). 

With  the  preceding  definitions,  the  power  spectrum  of  the  sam- 
pling error  for  two-dimensional  track  type  sampling  is  defined  as 

|Fg(u,v)|^  = |f^(u,v)|^  + 1f^(u,v)|^  (7) 


As  In  the  one-dimensional  case,  the  estimated  mean  square  sampling 
error  is  obtained  by  numerically  solving  the  two-dimensional  form  of 
Parseval*s  formula.  That  is,  for  an  NxN  grid  of  equally  spaced  values 
(Ax=Ay = one  data  interval)  of  a two-dimensional  function  f , the 
mean  square  sampling  error  is 


N=1  N=1 

— E E 
-.2  n=o  m=o 
N 


f^ 

n,m 


N/2-1  N72-1 


2 


(8) 


Figure  11  illustrates  the  relationship  which  exists  between  the 

two-dimensional  spatial  and  frequency  domains  for  track-type  surveys. 

Prafitical  application  of  this  track-sampling  theory  involves  the 

application  of  a two-dimensional  FFF,  with  appropriate  prewhitening, 
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and  the  numerical  solution  of  equations  (6)  - (8)  to  obtain  an  estimate 
of  the  mean  square  sampling  error  and  the  two-dimensional  sampling 
error  spectrum  at  any  desired  folding  frequency  which  is  a harmonic  of 
the  test  sample  field. 

As  a numerical  test  of  this  theory,  consider  the  contour  chart 

shown  in  Figure  12  to  be  the  result  of  a small  detailed  survey  within  a 

much  larger  homogeneous  province.  These  test  data,  with  units  of  arc 

seconds  (sec),  actually  consist  of  16x16  digital  values  derived  from  a 

model  of  the  north  component  of- vertical  deflection  at  a grid  spacing 

as  shown  in  Figure  12.  In  this  example,  the  y axis  is  assumed  to  be 

oriented  in  the  north  direction.  Figure  13  is  a contour  chart  of  all 

four  quadrants  of  the  two-dimensional  amplitude  spectrum  of  this  base 

data  after  appropriate  prewhitening  and  correction.  This  figure,  as 

well  as  subsequent  contour  charts  of  computed  two-dimensional  spectra, 

indicates  a region  centered  at  F = F = 0.0  in  which  the  spectral  con- 

X y 

tent  is  undefined.  This  is  a consequence  of  the  prewhitening  operation 
which  essentially  removes  wavelength  components  longer  than  the  data 
record,  and  the  result  of  applying  a two-dimensional  low-pass  smooth- 
ing filter  to  the  spectral  estimates  in  order  to  show  the  general 
character  of  the  spectrum. 

The  numerical  solution  of  equations  (6)  - (8)  for  north-south 
tracks  spaced  at  varying  multiples  of  the  original  data  spacing 
produced  estimates  of  the  sampling  error  as  shown  in  Figure  lA.  For 
compa’*i'5^n,  the  error  estimates  resulting  from  the  utilization  of  a 
spectrum  computed  without  prewhitening  are  also  shown.  In  this 
particular  example,  the  effect  of  spectral  leakage  on  the  estimated 
sampling  error  spectrum  is  not  particularly  pronounced. 
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Figure  12.  Base  data  for  two-dimensional  sampling  test — 16  X 16  point  grid. 
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Figure  13.  Two-dimensional  amplitude  spectrum  af  base  data. 
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Figure  14.  Numerical  two-dimensional  test — Sampling  error  vs.  N-S  track  spacing. 
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At  each  track  spacing,  the  solution  of  equations  (6)  and  (7)  pro- 
duce an  estimate  of  the  two-dimensional  spectrum  of  the  sampling  error. 


Figures  15-17  are  contour  charts  of  the  two-dimensional  amplitude 
spectrum  of  the  sampling  error  at  selected  track  spaclngs.  As  expected 
from  the  preceding  discussion  concerning  the  operations  of  S^(u,v),  and 
the  general  "red  noise"  type  character  of  the  spectrum  of  the  base 
data  (Figure  13),  Figures  15-17  indicate  that  the  maximum  amplitude  of 
the  sampling  error  spectrum  decreases  as  the  spacing  between  tracks 
decreases,  and  that  the  maximum  amplitude  is  located  quite  close  to  the 
folding  frequency. 

As  a numerical  test  of  the  accuracy  of  these  sampling  error  esti- 
mates, every  second  column  of  data  was  extracted  from  the  base  data 
grid  to  simulate  the  survey  data  which  would  be  collected  along  north- 
south  tracks  spaced  two  data  intervals  apart.  Values  of  the  field  at 
the  missing  grid  points  were  generated  by  the  application  of  a cubic 
spline  interpolation  procedure  (Davis  and  Kontis,  1970)  designed  specif- 
ically for  gridding  data  from  track  - type  surveys.  The  actual  mean 
square  error  computed  for  this  test  case  is  shown  in  Figure  14.  The 
two-dimensional  amplitude  spectrum  of  this  error  is  contoured  in 
Figure  18.  A comparison  of  Figures  16  and  18  indicates  the  generally 
excellent  agreement  between  the  theoretical  estimate  of  the  error 
spectrum  and  the  actual  spectrum  from  this  numerical  test.  Figure  19 
is  a profile  through  these  two  error  spectra  along  the  line 
The  slig>'*-  disagreement  bet^-^een  the  two  curves  at  the  lower  frequen- 
cies is  undoubtedly  generated  by  two-dimensional  spectral  leakage 
since  no  prewhitening  filter  was  applied  to  the  sampling  error 
produced  by  the  numerical  test. 
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Figure  1 5.  Two-dimensional  error  spectrum,  N-S  track  spocing=4.0  Dl. 
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Figure  16.  Two-dimensional  error  spectrum,  N-S  track  spocing=2.0  Dl. 
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Figure  17.  Two-dimensional  error  spectrum,  N-S  track  spacing=1.14  Dl. 
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Figure  1 8.  Actual  error  spectrum  for  o N-S  track  spacing  of  2.0  Dl. 
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Figure  19.  Profile  of  sampling  error  spectrum  along  Fx=Fy  from  numerical  test. 


43 


In  the  practical  application  of  this  theory,  these  base  data  are 


obtained  by  a detailed  survey  of  a small  test  area  located  within  a 
large  homogeneous  province.  Since  this  detailed  survey  should  be 
designed  to  adequately  sample  the  base  data  in  all  directions,  a simple 
rotation  of  the  geographic  direction  of  the  x-y  grid  will  allow  the 
survey  design  process  to  be  utilized  to  produce  contours  of  the  esti- 
mated sampling  error  as  a function  of  both  track  spacing  and  track 
direction. 

Effect  of  Pre-Filtering  on  Sampling  Error 

At  first  glance,  it  would  appear  reasonable  that  low-pass  filter- 
ing the  digital  data  collected  along  survey  tracks  would  be  an 
efficient  means  of  reducing  the  mean  square  sampling  error.  The  fact 
that  this  is  generally  not  the  case  may  be  shown  in  the  following 
manner . 

The  two-dimensional  filter  weight  function  W(x,y) , which  represents 
the  application  of  a one-dimensional  digital  filter  G(y)  to  the  track 
data  assumed  parallel  to  the  y axis,  is  given  by  W(x,y)  = <S(x)g(y). 

The  two-dimensional  Fourier  transform  of  W(x,y)  is 


W(u,v) 


5(x)g(y)e  ’■^“^'^^^dxdy 


G(v). 


Figure  20  is  a plot  of  this  model.  Notice  that,  although  in  practice 
the  sampling  function  S^(x,y)  is  applied  first  and  then  a one-dimensional 
filter  i'-  •’pplied  to  this  sampled  data,  this  is  identical  to  applying 
W(x,y)  first  and  then  sampling  the  result.  Thus,  in  this  case. 


Dp(x,y)  = (f  (x,y)*W(x,y))  S^(x,y)  = [ f (x,y) S,^(x,y) ]*U(x,y)  , (9) 
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Figure  20.  The  two-dimensional  equivalent  of  olong-trock  filtering. 
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where  D„(x,y)  is  the  filtered  data.  In  the  frequency  domain » the  right 
F 

hand  side  of  equation  (9)  is  equivalent  to  [F(u,v)*S^(u,v)]W(u,v) . 

Figure  21  illustrates  this  operation.  It  is  clear  that  the  net  result 
of  pre-filtering  the  track  data  is  to  introduce  an  error  of  omission 
for  frequencies  in  the  along-track  direction  while  having  no  effect  on 
either  the  error  of  commission  or  the  error  of  omission  in  the  cross 
track  direction.  For  frequencies  in  other  directions,  the  application 
of  W(u,v)  may  reduce  the  sampling  error  to  some  extent,  but  this  effect 
must  be  evaluated  on  a case  by  case  basis  in  order  to  determine  if  the 
net  reduction  in  the  error  of  commission  is  significant  in  light  of  the 
increase  in  the  error  of  omission. 

The  inherent  limitation  in  prefiltering  track  data  by  the  applica- 
tion of  a low-pass  filter  to  the  digital  data  should  not  be  construed 
to  mean  that  nothing  may  be  done  to  reduce  sampling  error,  On  the  con- 
trary, in  many  practical  situations  it  is  possible  to  exert  a consider- 
able amount  of  control  over  this  error.  This  control  may  be  accomplish- 
ed either  through  the  use  of  the  inherent  physical  properties  of  the 
data  which  are  being  sampled  or  through  the  control  of  the  two- 
dimensional  frequency  response  of  the  measuring  device.  For  example, 
in  airborne  surveying  of  magnetic  or  gravity  fields,  the  natural 
exponential  decay  of  the  two-dimensional  spectrum  as  a function  of 
height  above  the  source  (Grant  and  West,  1965)  results  in  the  natural 
application  of  a two-dimensional  low-pass  filter.  As  the  flight 
elevation  is  increased,  the  amplitude  of  the  short  v.’avelength  compon- 
ents t-he  field  are  reduced.  Thus,  while  the  error  of  omission  is 
increased,  the  error  of  commission  v</ill  be  appreciably  reduced.  This 
results  in  a relatively  undistorted  measurement  of  the  longer  wave- 
length components  of  the  field. 
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Figure  21.  The  effect  of  olong-trock  filtering  on  sampling  error. 
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An  example  of  the  control  of  sampling  error  through  the  design  of 
instrumentation  may  be  found  in  the  measurement  of  gcoid  height  by  the 
use  of  a satellite  mounted  radar  altimeter.  Adjusting  the  beam  width 
of  the  receiver  controls  the  size  of  the  illuminated  spot  on  the  sea 
surface  resulting  in  a two-dimensional  filtering  operation  which 
reduces  the  error  of  commission  at  the  expense  of  increasing  the  error 
of  omission. 
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ALGORITHMS  FOR  DETERMINING  HOMOGENEOUS  PROVINCES 


The  successful  practical  application  of  the  theory  of  survey  design 
depends  to  a great  extent  upon  the  accuracy  with  which  the  determination 
of  the  boundaries  of  relatively  homogeneous  provinces  may  be  made.  In 
a statistical  sense,  a stochastic  process  is  defined  as  being  stationary, 
in  the  wide  sense,  if  its  mean  and  autocorrelation  function  are  indepen- 
dent of  position  (Papoulis,  1962).  This  is  equivalent  to  requiring  that 
the  two-dimensional  amplitude  or  power  spectrum  be  independent  of  posi- 
tion. It  is  clear  from  the  preceeding  description  of  the  way  in  which 
the  track  sampling  model  (S^(u,v))  operates  on  the  Fourier  transform  of 
the  data,  that  the  requirement  for  stationarity  may  be  somewhat  relaxed 
for  survey  design.  Within  this  context,  a province  is  defined  as  homo- 
geneous if  the  sample  data  are  nearly  stationary  only  for  those 
frequency  components  that  will  contribute  to  sampling  error  at  a partic- 
ular track  spacing. 

Even  with  this  relaxation  of  the  requirement  for  stationarity,  a 

rigid  adherence  to  this  definition  is  not  practically  feasible  since 

this  would  require  that  different  province  boundaries  be  defined  for 

each  track  spacing  and  each  track  direction.  The  particular  process 

which  is  utilized  to  define  these  province  boundaries  depends  upon 

several  factors  such  as  the  availability  and  extent  of  reconnaissance 

data  and  whether  the  data  is  basically  a function  of  two-dimensions 

such  as  gravity  and  bathymetry,  or  three  dimensions  such  as  oceanic 

temperature  and  salinity.  In  addition,  the  appropriate  orientation  of 

the  survey  lines  may,  in  some  instances,  be  controlled  by  factors  other 

than  the  local  spectral  content  which  is  normally  used  to  delineate 

provinces.  One  example  of  this  is  the  effect  of  errors  in  the 
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determination  of  heading  and  velocity  for  computing  the  Eotvos  correc- 
tion to  be  applied  to  shipboard  gravity  survey  data.  In  this  case,  the 
Eotvos  correction  contains  less  error  if  the  survey  tracks  are  oriented 
east-west  (Glicken,  1962) • The  following  discussion  consists  of  a 
detailed  description  of  the  numerical  algorithms  which  have  been  devel- 
oped for  determining  the  location  of  the  boundaries  of  homogeneous 
provinces  through  the  use  of  reconnaissance  data.  Fortran  IV  computer 
programs  are  available  from  the  author  for  implementing  each  of  these 
algorithms. 


Two-Dimensional  Fields 

In  general,  the  reconnaissance  data  which  are  available  for  delin- 
eating provinces  for  gravity,  magnetic,  or  bathymetric  surveys  consist 
of  widely  spaced  tracks  of  generally  random  orientation.  Since  the 
application  of  S^(u,v)  results  in  a replication  of  the  transform  of  the 
data  in  the  cross-track  direction,  it  is  desirable  that  a sufficient 
amount  of  reconnaissance  data,  or  related  geophysical  information,  be 
available  for  determining  the  general  direction  of  any  major  lineation 
of  the  field.  In  order  to  minimize  the  mean  square  sampling  error,  the 
most  efficient  orientation  of  the  final  survey  tracks  is  in  a direction 
normal  to  these  lineations.  This  assumes,  of  course,  that  there  is  no 
overriding  constraint  on  track  direction  such  as  that  imposed  by  the 
Eotvos  correction.  The  object  of  this  province  selection  process  is  to 
outline  regions  in  which  short  wavelength  energy  in  the  cross-track 
directio-'  is  expected  to  be  uniform  based  on  the  reconnaissance 
information. 

At  the  present  time,  the  algorithm  which  has  been  developed  for 
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this  type  of  province  selection  consists  of  the  following  steps: 

1.  Those  reconnaissance  tracks  which  arc  most  nearly  parallel 
to  the  major  lineations  are  selected  for  analysis.  Each 

of  these  tracks  is  separated  into  relatively  straight  line 
segments  and,  following  the  procedure  of  Davis  and  Kontis 
(1970) , a least-squares  straight  line  is  fitted  to  the 
positions  of  the  data  points  in  each  segment.  The  actual 
position  of  each  data  value  is  then  mapped  onto  this 
straight  line  and  the  data  value  is  adjusted  by  utilizing 
the  local  gradient  of  the  field.  A cubic  spline  is  then 
applied  to  interpolate  new  data  values  which  are  equally 
spaced  in  terms  of  dis'cance  down-track.  The  numerical 
details,  as  well  as  tests  indicating  the  accuracy  of  this 
procedure,  are  given  in  Davis  and  Kontis  (1970). 

2.  At  this  point,  a practical  decision  must  be  made  concerning 
the  maximum  track  spacing  which  one  would  be  willing  to 
allow  in  the  final  survey  design.  This  decision  is  based 
on  an  estimate  of  sampling  error  as  a function  of  sample 
spacing  obtained  by  applying  equations  (2)  and  (3)  to  a 
segment  of  reconnaissance  data  selected  in  step  (1).  This 
segment  of  data  should  be  selected  to  contain  a relatively 
small  amount  of  short  wavelength  energy.  This  maximum 
track  spacing  defines  a folding  frequency  which  is  used  in 
conjunction  with  equations  (A-1)  and  (A-2) , shown  in 
Appendix  A,  to  design  a high-pass,  one-dimensional,  digital 
filter  which  will  pass  only  those  frequencies  that  will 

contribute  to  sampling  error. 
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3.  The  equally  spaced  reconnaissance  data  produced  in  step  (1) 
are  high-pass  filtered  with  the  filter  designed  in  step  (2), 
and  an  approximation  to  the  envelope  is  computed.  The 
variability  of  the  amplitude  of  this  envelope  is  used  as 
a practical  means  for  locating  the  boundaries  of  stationary 
data  segments. 

Although  the  means  for  computing  a mathematically 
correct  envelope  via  the  frequency  domain  application  of  a 
Hilbert  transform  (Papoulis,  1962)  are  available  on  a large 
scale  computer,  an  approximation  is  used  here  for  applica- 
tion to  small  shipboard  computers.  This  approximation  con- 
sists simply  of  computing  the  absolute  value  of  the  high- 
passed  data  and  applying  a low-pass  filter  with  a very  low 
frequency  cutoff  and  a small  number  of  v;eights  to  conserve 
data. 

A.  In  order  to  generate  the  locations  of  the  province  boundaries 
from  the  envelope  estimates,  two  additional  decisions  must 
be  made  at  this  point.  These  are: 

A.  What  is  the  minimum  size,  ie.  horizontal  extent,  of  a 
province  which  one  would  be  willing  to  accept  in  terms  of 
modifying  the  track  spacing? 

B.  In  terras  of  the  envelope  of  the  high-passed  data,  how 
much  variability  will  be  allowed  within  a province? 

While  these  are  obviously  rather  arbitrary  decisions,  in 

general,  the  appropriate  values  of  these  control  parameters 

usually  become  apparent  after  several  test  runs  on  the 

reconnaissance  data.  As  a point  of  reference,  the  present 
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minimum  province  size  selected  for  shipboard  gravity 
survey  design  is  32  nautical  miles  (WI)  with  an  allow- 
able envelope  variability  of  1.5  mgals. 

5.  Utilizing  the  above  control  parameters,  provinces  are 
selected  which  are  of  at  least  minimum  size  and  within 
which,  the  envelope  variability  is  within  the  acceptable 
range.  If  reconnaissance  tracks  are  available  at 
various  orientations,  these  steps  may  be  repeated  to 
delineate  provinces  for  any  desired  orientation  of  the 
final  survey. 

Given  these  province  boundaries  estimated  from  the  selected  recon- 
naissance data  segments,  the  final  steps  in  the  survey  design  procedure 
consist  of  conducting  a detailed,  closely-spaced,  track-type  survey  of 
a small  test  area  within  each  province  and  applying  the  previously 
defined  procedure  to  obtain  the  design  parameters  for  each  type  province. 
These  design  parameters,  or  decision  products,  consist  of  the  mean 
square  sampling  error  as  a function  of  track  spacing,  as  well  as  track 
direction  if  appropriate,  and  the  two-dimensional  spectral  content  of 
the  sampling  error  at  selected  track  spacings.  Combining  these  error 
estimates  with  a predetermined  accuracy  specification  results  in  the 
final  survey  design.  As  a consequence  of  having  available  the  spectral 
content  of  the  sampling  error,  it  is  possible  to  accommodate  a wide 
range  of  accuracy  specifications.  For  example,  if  an  accuracy  specifi- 
cation is  stated  in  terms  of  the  allowable  mean  square  sampling  error 
within  a particular  frequency  band,  the  application  of  equation  (8)  to 
the  sampling  error  spectrum  over  this  restricted  frequency  band  will 
produce  the  required  estimates  of  sampling  error.  Figure  22  is  a 
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RECONNAISSANCE  SURVEY  (WIDE  TRACK  SPACING) 


GENERAL  PROVINCE  DESIGNATION  AND  SELECTION  OP  TYPE  AREAS 


DECISION  PRODUCTS  FOR  EACH  TYPE  AREA 
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1. 

Track 
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2. 
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3. 

Down-Track  Scunple  Rate 

Figure  22.  Survey  design  procedures  for  two-dimensional  fields. 
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general  flow  diagram  outlining  the  procedure  which  is  utilized  for 
designing  oceanographic  surveys  of  gravity,  magnetics,  or  bathymetry. 

One  additional  point  should  be  made  with  regard  to  selecting  the 
proper  track  spacing  for  the  detailed  type  area  surveys.  At  the 
present  time,  the  most  appropriate  procedure  appears  to  consist  of 
selecting  a section  of  reconnaissance  track  data  which  exhibits  the 
greatest  magnitude  of  short  wavelength  energy  within  each  selected 
province.  A two-dimensional  Fourier  transform  or  equivalently,  a 
Hankel  transform,  is  applied  to  this  selected  data  set  which  is  assumed 
to  possess  circular  symmetry.  This  transform  is  then  used  with  equa- 
tions (7)  and  (8)  to  select  a track  spacing  for  the  type  area  survey 
such  that  a minimum  amount  of  sampling  error  will  be  generated  in  the 
type  area  data. 

Three-Dimensional  Fields 

The  process  x-?hich  has  been  developed  for  determining  the  boundaries 
of  homogeneous  provinces  for  three-dimensional  fields  such  as  oceanic 
temperature,  salinity,  or  sound  speed  is  considerably  more  complex  than 
the  two-dimensional  case.  It  should  be  pointed  out  that  these  fields 
are,  in  reality,  four-dimensional,  since,  as  will  be  shown  in  a later 
example,  temporal  variability  is  quite  pronounced  in  some  ocean  areas. 
Despite  this  fact,  the  present  lack  of  appropriate  data  restricts  the 
survey  design  procedure  to  consider  only  three  dimensions. 

As  the  basis  for  this  process,  a generalized  procedure  has  been 
developed  for  modeling  three-dimensional  fields  through  the  utiliza- 
tion of  reconnaissance  data  consisting  of  randomly  spaced  vertical  pro- 
files. The  historical  data  which  are  available  for  models  of  oceanic 
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temperature,  salinity,  and  sound  speed  are  processed  in  the  form  of  a 
digital  data  bank  by  the  National  Oceanographic  Data  Center  (NODC, 

196A).  The  pertinent  infonnation  in  this  data  bank  consists  of  the 
latitude,  longitude,  and  time  of  each  ocean  station  along  with  the 
recorded  values  of  temperature  and  salinity,  and  the  computed  values  of 
sound  speed  at  what  are  defined  as  standard  oceanographic  depths.  As  a 
consequence  of  the  large  vertical  gradients  of  these  fields  at  the 
shallower  depths,  the  vertical  spacing  between  samples  normally  varies 
as  a function  of  depth.  These  standard  oceanographic  depths  were  deter- 
mined by  the  International  Association  of  Physical  Oceanography  in  1936 
(Sverdrup  et  al,  19A2)  and  are  presented  in  Table  (1). 

Table  1.  Standard  Oceanographic  Depths 


Depth  (n) 


10 

75 

300 

800 

2000 

20 

100 

AOO 

1000 

2500 

30 

150 

500 

1200 

3000 

50 

200 

600 

1500 

and  every 

1000 

thereafter 


At  each  standard  depth,  sound  speed  is  computed  via  the  empirically 

derivf’d  '-'’lationship  equating  sound  speed  to  the  measured  values  of 

temperature,  salinity,  and  pressure  (h’ilson,  1960). 

Taking  sound  speed  as  an  example,  a primary  three-dimensional  model 

is  constructed  in  the  following  manner.  To  reduce  the  effect  of 
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seasonal  variations,  the  station  data  are  separated  into  appropriately 
defined  winter  and  summer  stations.  V.’orking  with  one  vertical  profile 
at  a time,  a cubic  spline  interpolation  procedure  is  applied  to  the 
sound  speed  values  at  the  standard  depths  to  generate  data  values  at  an 
equal  depth  spacing  of  50  meters.  A seven  degree  least  squares  orthog- 
onal polynomial  (Van  Voorhis  and  Davis,  1964)  is  then  fitted  to  these 
equally  spaced  profile  values  over  the  depth  range  200  to  2450  meters. 
The  starting  depth  of  200  meters  was  selected  because  of  the  high 
frequency  temporal  variations  which  normally  occur  at  the  shallower 
depths.  Extensive  tests  indicate  that  a seven  degree  polynomial  will 
generally  yield  an  RMS  residual  of  less  than  1 m/sec.  This  level  of 
residual  was  selected  on  the  basis  of  a study  by  Pickett  (1972)  which 
indicates  that  the  error  in  computed  sound  speed  is  within  0.8  to  2.8 
m/sec  depending  on  salinity  and  depth. 

Taking  the  first  coefficient  of  the  polynomial  thus  determined  for 
each  profile  within  a selected  area,  a two-dimensional  gridding  algo- 
rithm is  applied  to  interpolate  the  value  of  this  coefficient  at 
equally  spaced  grid  points.  As  a consequence  of  the  paucity  of  existing 
data,  the  grid  spacing  for  sound  speed  is  presently  selected  to  be  30 
minutes  of  latitude  x 30  minutes  of  longitude.  This  two-dimensional 
gridding  algorithm,  which  was  developed  by  Rankin  (1974),  is  essential- 
ly a three  phase  algorithm  consisting  of  local  averaging,  distance 
weighting  (Shepard,  1963),  and  cubic  spline  interpolation.  Applying 
this  interpolation  procedure  to  each  of  the  eight  polynomial  coeffi- 
cients derived  at  each  station  produces  a set  of  eight  coefficient 
surfaces  which  comprise  the  primary  three-dimensional  sound  speed  model. 

In  order  to  develop  province  boundaries  from,  this  numerical  model, 
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a two-dimensional  low-pass  filter  is  applied  to  the  grlddod  coefficient 


surfaces  for  degree  one  through  four.  The  weight  function  for  this 

filter  is  derived  from  equation  (A-3)  of  Appendix  A with  a K of  100 

c 

data  intervals  and  a of  5 data  intervals.  For  a 30x30  minute  grid 
spacing,  this  produces  a filter  which  effectively  removes  wavelengths 
shorter  than  2.5  degrees  of  latitude  or  longitude,  thus  reducing  the 
effect  of  steep  local  horizontal  gradients  in  the  coefficient  surfaces. 
These  local  gradients  are  normally  generated  by  short  term  temporal 
variations  occurring  at  closely  spaced  stations,  or  by  erroneous  data. 

By  restricting  the  two-dimensional  frequency  content  of  the  polynomial 
surfaces  in  this  way,  the  province  boundaries  are  defined  on  more  of  a 
regional  basis  than  on  a local  basis.  Only  the  degree  one  through  four 
coefficient  surfaces  are  used  to  define  the  sound  speed  province  bound- 
aries since  these  have  been  found  to  account  for  most  of  the  variability 
in  acoustic  propagation  over  this  depth  range. 

The  next  step  in  the  process  consists  of  applying  a two-dimensional 
bicubic  spline  (Davis  and  Kontis,  1970)  to  each  of  the  filtered 
coefficient  surfaces  in  order  to  determine  the  partial  derivatives 
(•|^,  > 3nd  the  magnitude  of  the  gradient  at  each  grid  point.  As  was 

the  case  for  provinces  defined  for  two-dimensional  fields,  a decision 
must  be  made  as  to  the  minimum  acceptable  size  of  a province.  Uti- 
lizing this  minimum  province  size,  an  appropriate  magnitude  of  the 
gradient  is  selected  from  the  contoured  gradient  values  to  define  the 
province  boundaries  for  each  coefficient  surface.  A composite  of  these 
four  boundary  charts  produces  the  final  estimates  of  what  are  defined 
as  secondary  province  boundaries  for  survey  design.  Figure  23  is  a 

general  flow  diagram  for  this  province  selection  algorithm. 
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Figure  23.  Algorithm  for  defining  sound  speed  province  boundaries. 
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In  sununary,  the  present  technique  for  selecting  province  boundaries 
for  the  three-dimensional  oceanic  fields  is  based  on  locating  thbse 
areas  in  which  the  overall  character  of  a vertical  profile  changes 
abruptly  in  the  horizontal  direction.  This  is  in  contrast  to  the  two- 
dimensional  case  in  which  the  province  boundaries  are  defined  by  a 
prescribed  change  in  the  magnitude  of  the  high  frequency  spectral 
content. 
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PRACTICAL  APPLICATIONS 


Several  practical  applications  are  presented  in  order  to  demon- 
strate the  adaptability  of  the  preceeding  theory  of  survey  design,  and 
the  efficiency  with  which  the  two-dimensional  frequency  domain 
representation  of  sampling  error  may  be  propagated  through  various 
linear  operations. 

The  application  of  survey  design  to  hydrographic  surveying  illus- 
trates the  utilization  of  a one-dimensional  modification  of  the  two- 
dimensional  theory  for  use  with  small  scale  computers  in  near 
real-time.  This  modification  is  possible  if  the  appropriate  orienta- 
tion of  the  survey  tracks  is  known  a priori,  and  the  desired  accuracy 
specification  is  in  terms  of  a broad-band  mean  square  or  RMS  error 
rather  than  the  spectral  content  of  the  error. 

The  application  of  the  two-dimensional  theory  to  the  design  of 
gravity  surveys  for  computing  vertical  deflection  and  geoid  undulation 
was  selected  to  illustrate  one  of  the  more  intriguing  and  useful 
aspects  of  the  theory.  In  this  application,  the  desired  accuracy  speci- 
fication is  stated  in  terms  of  both  a broad  band  mean  square  error,  and 
the  two-dimensional  spectral  content  of  the  error  in  geoid  undulation 
as  well  as  the  error  in  the  north  and  east  components  of  vertical 
deflection.  The  procedure  for  satisfying  this  type  of  specification 
consists  of  propagating  the  two-dimensional  spectrum  of  the  sampling 
error  generated  by  the  basic  gravity  survey  through  the  linear 
operations  required  for  computing  undulation  and  deflection  components. 
This  application  of  the  survey  design  theory  illustrates  the  accuracy 
and  inherent  efficiency  of  this  frequency  domain  procedure  as  well  as 

the  effect  of  regional  lineations  on  sampling  error. 
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The  application  of  the  two-diraensional  survey  design  theory  to 
three-dimensional  fields  is  illustrated  by  an  example  of  a survey 
design  for  oceanic  sound  speed.  This  application  also  Illustrates  the 
facility  with  which  this  sampling  theory  may  be  adapted  to  provide  a 
very  useful  estimate  of  the  sampling  error  as  a function  of  depth. 

One  final  practical  application  of  the  two-dimensional  theory 
Illustrates  the  effect  of  sampling  error  on  a classical  problem  of 
geophysical  interpretation.  This  example  was  specifically  designed  to 
show  that  an  apparently  insignificant  amount  of  sampling  error  can 
severely  influence  an  Interpretation  if  the  interpretation  technique  is 
sensitive  to  the  high-frequency  components  of  the  anomalous  field. 

Near  Real-Time  Hydrographic  Survey  Design 

At  the  present  time,  the  U.  S.  Naval  Oceanographic  Office  is  in 
the  process  of  developing  a high-speed  hydrographic  survey  and  charting 
system  (HYSURCH)  for  defining  the  topography  of  the  ocean  bottom  in 
shallow  coastal  areas.  This  automated  system  is  designed  to  enable  the 
Navy  to  rapidly  and  economically  conduct  hydrographic  surveys  to  a pre- 
determined accuracy  as  well  as  on-site  data  evaluation  and  automated 
field  production  of  contour  charts.  The  basis  for  the  economy  of  this 
approach  is  the  near  real-time  application  of  survey  design  procedures 
in  the  following  manner. 

Figure  24  shows  the  general  concept  which  is  essentially  a one- 
dimensional modification  of  the  two-dimensional  sampling  theory 
developed  for  the  design  of  gravity,  magnetic,  and  deep  ocean  bathy- 
metric surveys.  The  major  regiona'l  trends  of  the  field  normally  are 

parallel  to  the  coastline.  On  this  basis,  the  appropriate  direction  of 
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RECONNAISSANCE  TRACKS  PARALLEL  TO  REGIONAL  TREND 
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ONE-DIMENSIONAL  ERROR  ANALYSIS  AND  SURVEY  DESIGN 


Figure  24.  Near  real-time  hydrographic  survey  design  (HYSURCH). 


the  survey  tracks  is  normal  to  the  coastline.  With  this  constraint  on 
the  survey  orientation,  data  are  collected  along  reconnaissance  lines 
run  generally  parallel  to  the  coastline.  The  digital  data  collected 
along  these  lines  is  then  used  to  define  the  locations  of  homogeneous 
province  boundaries  utilizing  the  algorithm  developed  for  application 
to  two-dimensional  fields. 

A one-dimensional  Fourier  transform  is  applied  to  the  data  within 
each  of  these  selected  provinces  for  each  of  the  available  reconnais- 
sance tracks.  With  these  spectral  estimates,  equations  (3)  and  (A)  aye 
numerically  evaluated  to  obtain  an  estimate  of  the  sampling  error  as  a 
function  of  sample  spacing.  Since  the  actual  survey  tracks  will  be 
oriented  normal  to  these  reconnaissance  lines,  this  sampling  error  is, 
in  reality,  a function  of  the  survey  track  spacing.  Given  a required 
survey  accuracy  specification,  the  appropriate  sample/ track  spacing  is 
selected  for  each  defined  province  and  the  final  survey  is  completed 
using  the  average  specified  track  spacing  for  each  area. 

At  the  conclusion  of  the  survey  operation,  the  track  data  are 
examined  to  determine  if  any  significant  provinces  v;ere  missed  by  the 
preliminary  reconnaissance  tracks.  In  this  case,  additional  short 
survey  lines  will  be  run  as  necessary  and  the  sampling  error  will  be 
estimated  at  the  actual  track  spacing  achieved  within  each  of  these 
additional  provinces. 

Figure  25  is  an  example  of  a test  of  the  province  selection  process 
applied  to  survey  data  collected  along  one  reconnaissance  track.  The 
lower  curve  is  a plot  of  the  original  sounding  data  collected  at  five- 
meter  intervals  along  the  track.  The  upper  curve  is  a plot  of  the 
original  data  after  high-pass  filtering  with  a frequency  response  set  to 
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Figure  25.  Example  of  HYSURCH  province  selection  process. 


pass  spectral  components  which  would  contribute  to  the  sampling  error 
at  a maximum  allowable  track  spacing  of  2000  meters.  A plot  of  the 
approximation  to  the  envelope  of  this  high  frequency  energy  is  also 
shown.  Selecting  a minimum  province  size  of  160  meters  in  the  hori- 
zontal dimension,  and  a maximum  variability  of  the  envelope  within  a 
province  of  0.2  meters,  the  province  selection  algorithm  produced  posi- 
tions for  the  province  boundaries  as  shown  by  the  dashed  lines.  On 
this  particular  test,  there  were  three  type-one  provinces  selected  with 
envelope  amplitudes  within  the  range  0-0.2  meters,  and  two  type-two 
provinces  with  amplitudes  in  the  range  0.2-0. 4 meters.  Figure  26  is  a 
plot  of  the  estimated  RMS  of  the  sampling  error  as  a function  of 
sample/track  spacing  for  the  first  province  type  one  segment.  The  geo- 
graphic coordinates,  in  x-y  units,  of  this  province  are  also  noted. 
Figure  27  shows  the  relationship  which  exists,  for  this  data  set, 
between  the  FFT  amplitude  spectrum  computed  on  the  raw  data  without 
window  modification  or  prewhitening,  the  spectrum  of  the  prewhitened 
data,  the  rough  spectrum  corrected  for  the  prewhitening  operation,  and 
the  final  smoothed  corrected  spectrum.  As  was  the  case  for  the  gravity 
example  shown  in  Figure  3,  the  requirement  for  prewhitening  is  readily 
apparent. 


Survey  Design  for  Vertical  Deflection 
and  Geoid  Undulation 

In  this  example,  the  requirement  is  to  design  a track-type  survey 
operation  to  produce  products  of  a specified  accuracy  in  which  the 
accuracy  requirement  is  placed  not  on  the  data  actually  being  collect- 
ed, but  upon  two  quantities  (vertical  deflection  and  geoid  undulation) 
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Figure  26.  Sompling  error  vs.  sample/track  spacing  for  first  province  segment. 
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FFT  AA^ITUOE  SPECTRUM  (METERS  * Dl ) 


Figure  27.  Amplitude  spectrum  of  data  from  first  province  segment. 
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computed  from  these  data.  In  addition,  although  the  overall  broad-band 
mean  square  sampling  error  is  of  general  interest,  particular  emphasis 
is  placed  on  estimating  the  two-dimensional  spectral  content  of  the 
sampling  error  in  order  to  provide  a means  for  determining  the  mean 
square  error  of  the  computed  quantities  within  a specified  frequency 
band. 

At  the  present  time,  the  vast  majority  of  existing  values  of  verti- 
cal deflection  and  geoid  undulation  are  computed  utilizing  worldwide 
measurements  of  the  earth's  gravity  field  and  the  mathematical  techni- 
ques of  physical  geodesy.  Deflection  of  the  vertical  is  defined,  at  a 
point,  as  the  angular  difference  between  the  normal  to  an  equipotential 
surface  of  the  earth's  gravity  field,  and  the  normal  to  a smooth  mathe- 
matically defined  reference  ellipsoid  having  the  same  potential 
(Heiskanen  and  Moritz,  1967).  This  vertical  deflection  is  normally 
computed  in  teirms  of  a north-south  component  (C)  and  an  east-west  com- 
ponent (n) . The  geoid  undulation  is  defined  as  the  distance  between 
these  two  equipotential  surfaces  measured  along  a line  normal  to  the 
reference  ellipsoid. 

In  order  to  compute  geoid  undulation  and  the  components  of  vertical 
deflection  from  measurements  of  the  gravity  field,  Stokes  integral 
(Heiskanen  and  Moritz,  1967)  is  numerically  solved  to  produce  the  undul- 
ation estimate,  and  the  Vening  Melnesz  integral  (Vening  Meinesz,  1928) 

Is  utilized  to  produce  the  values  of  the  vertical  deflection  components. 
In  each  of  these  integrals,  the  integrand  is  infinite  at  the  origin  of 
the  computation  system  which  is  centered  over  the  point  at  which  the 
deflection  or  undulation  estimate  is  desired.  This  presents  some  numer- 
ical difficulties  which  require  special  techniques  in  order  to  evaluate 
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this  so-called  Inner-zone  contribution  (Kontis,  et  al,  197^).  Excluding 
this  inner-zone  of  some  radius  the  Vening  Meinesz  integral  in 
spherical  polar  coordinates  (R,i|;,a)  (Figure  28)  is 
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where  and  Hp  are  the  north  and  east  deflection  components,  in  sec, 
at  the  computation  point,  p,  y is  the  mean  value  of  the  earth's  gravity 
field,  is  the  angular  distance  between  the  radius  vector  to  p and  a 
surface  element  of  area  R^  sin  i>)dii»da,  R is  the  mean  radius  of  the  earth, 
and  a is  the  azimuth,  measured  from  north,  of  the  surface  element  relative 

to  p.  Ag(^i),a)  is  the  mean  value  of  free-air  gravity  in  each  surface  ele- 
ment, and  is  the  Vening  Meinesz  function  given  by 


V(ip) 


. -cos  ij;/2 
_2sin^i|;/2 


- 6 cos  )\)/2  + 8 sim{>  +3  sin\|^  ln(sin^i|//2+sini{>/2) 


-3 (1-sin  (i>>/2) 
sinip 


Utilizing  the  same  system  of  spherical  polar  coordinates,  and  exclud- 
ing the  inner-zone,  the  Stokes  integral  for  computing  undulation  (Np)  is 


N = 


p 2ttyR 


■'a=o 


Ag(i^,a)S(\i;)R^  sin  ii;di|jda. 


(11) 


where  S(ti<)  is  the  Stokes  function  given  by 


S(<^)  = ^ 


sin^t;j/2  ~ ^ 2^  ^ cos;^ln(sin^  ilj/2+sin  i{//2) 


Since  the  magnitude  of  both  V(i|;)  and  S(ii/)  decreases  rapidly  as 
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Figure  28.  Spherical  polar  coordinate  system. 
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Increases,  the  gravity  field  close  to  the  computation  point  has  a much 
larger  effect  on  the  computed  undulation  and  deflection  than  does  the 
gravity  field  at  greater  distances.  As  a result,  as  the  computation 
point  (p)  is  moved  about  to  obtain  a closely  spaced  grid  of  computed 
values,  the  high  frequency  spectral  components  of  deflection  and  undula- 
tion are  generated  primarily  by  the  local  short  wavelength  components 
of  the  gravity  field.  As  a consequence,  both  and  S(ip)  may  be 

approximated  locally  by  the  so-called  flat-earth  approximation 
(Heiskanen  and  Moritz,  1967).  This  approximation  for  VCij!)  may  be  deriv- 
ed as  follows.  On  a sphere  of  radius  R,  the  arc  distance  between  any 
two  points  is  R For  large  R and  small  i|;,  the  arc  distance  may  be 
approximated  by  the  linear  disr'ance  r and  the  spherical  surface  may  be 
approximated  by  a plane.  For  small  if),  the  first  term  in  V(ij))  is  dominant 
so  that 


With  the  element  of  surface  area  on  the  plane  given  by  rdrda,  the  flat- 
earth  approximation  to  equation  (10)  is  given  by 


A sin^if)/2  if)^ 


- CSC  1” 

2ity 


rdrda. 


In  cartesian  coordinates,  with  y directed  north-south, 
y X 

rdrdo  = dxdy,  cosa  = sin  a = — , and 

r r 


(12) 


where  r = /x?4y2  . 
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In  a similar  manner,  (Heiskanen  and  Moritz,  1967)  a flat  earth 
approximation  may  be  derived  for  gcoid  undulation.  This  approximation, 
in  polar  coordinates,  is  given  by 

~ 2^  I I cartesian 

coordinates,  by  N ~ j j dxdy  (13) 

It  is  apparent  that,  as  the  location  of  the  computation  point  (p)  is 

varied,  both  equations  (12)  and  (13)  describe  a two-dimensional  cross- 

correlation  process  between  the  free-air  gravity  observations  and  the 

y X X 

two-dimensional  operators  and  — . This  observation  has  been 

utilized  as  the  basis  for  a very  efficient  algorithm  for  computing 
deflection  and  undulation  in  a large  scale  production  operation 
(Mlchlik,  1973). 

This  process  is  more  readily  applicable  to  the  survey  design 

problem  when  the  cross-correlation  is  converted  into  a two-dimensional 

convolution.  This  is  simply  accomplished  by  observing  that  the  C 

operator  (^)  and  the  n operator  (^)  are  both  even  functions  with 

respect  to  one  coordinate  axis  and  odd  with  respect  to  the  other,  and 

that  the  undulation  operator  (•^)  possesses  circular  symmetry.  Thus, 

equation  (12)  may  be  changed  to  a convolution  operation  by  simply 

changing  the  sign  and  equation  (13)  requires  no  modification.  In  this 

form,  the  convolution  theorem  (Bracewell,  I965)  may  be  used  with  the 

Fourier  transform  of  these  operators  to  convert  the  two-dimensional 

sampling  error  spectrum  of  the  gravity  survey  data  into  the  sampling 

error  spectra  for  undulation  and  the  vertical  deflection  components. 

Application  of  this  process,  along  with  the  solution  of  the 
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two-dimensional  form  of  Parscval's  relationship,  results  in  a straight- 
forward design  of  a gravity  survey  which  will  produce  a specified 
error  in  the  geoid  undulation  and  the  components  of  vertical  deflection. 

The  two-dimensional  Fourier  transform  of  the  deflection  operators 
have  been  developed  by  Shaw  et  al  (1969)  and  are  given  by 


i CSC  1" 


(f  Vf 

- X y 


(lA) 


and 


n (f  »f  ) = 

o X y 


i CSC  1" 


L(f  ^+f 

^ X y 


(15) 


A contour  of  the  first  quadrant  of  n (f  ,f  ) is  shown  in  Figure  29.  A 

o X y 

90  degree  rotation  of  the  axes  in  Figure  29  produces  a contour  of  the 

second  quadrant  of  5^(f^,f^).  Since  the  geoid  undulation  operator 

possesses  circular  symmetry,  the  two-dimensional  Fourier  transform 

reduces  to  a Hankel  transform  (Bracewell,  1965).  Thus,  with 

N (x,y)  = — — , then 
o ^ 2TTYr 


N_(u,v)  = 2u 

o 


rJ 

I r 


[r(uW)^^^l 


O''  J ^ 2iTYr 

r(u^+v^)'^'  ‘ Idr 


Y(u*+v^) 


(Gradshteyn  and  Ryzhik,  1965) 


or,  converting  to  normalized  frequency. 


N (f  ,f  ) = 
o X y 


K 


2TiY(f  ^+f  ^)^^^ 

' X y 


(16) 
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Figur*  29.  First  quadrant  of  two-dim«nsior>al  frtquanqr  raspons*  of  tfw  ETA  oporotor: 


75 


where  K is  the  length  of  a data  interval  in  meters. 

As  an  illustration  of  the  accuracy  of  the  flat-earth  approximations 
and,  consequently,  the  accuracy  of  the  spectral  domain  conversion  of  the 
sampling  error  in  the  gravity  field  to  error  in  undulation  and  deflec- 
tion components,  a theoretical  gravity  reference  surface  was  developed 
from  a distribution  of  16  point  masses.  The  formulation  for  this 
process  is  given  in  Davis  and  Kontis,  1970.  Figure  30  is  a contour 
chart  of  this  reference  gravity  field  at  a contour  interval  of  5 mgals. 
The  gravity  field  was  computed  on  a one  minute  grid  covering  the  entire 
80x80  minute  area.  In  addition,  the  magnitude  and  direction  of  the 
anomalous  gravity  vector  was  utilized  to  compute  the  ? and  q deflection 
components  at  each  grid  point.  Bruns  formula  (Heiskanen  and  Moritz, 
1967),  together  with  the  values  for  the  anomalous  potential,  was  used 
to  determine  the  undulation  reference  surface. 

Figure  31  is  a contour  chart  of  the  appropriately  prewhitened  and 
corrected  two-dimensional  amplitude  spectrum  of  this  gravity  model 
field.  Figures  32-3A  are  contour  charts  of  the  actual  two-dimensional 
amplitude  spectra  for  the  gridded  C.  D and  N reference  model  surfaces. 
Figures  35-37  were  generated  by  applying  the  frequency  responses  of  the 
flat-earth  operators  directly  to  the  two-dimensional  spectrum  of  the 
reference  gravity  field.  Except  for  the  high  frequency  roundoff  error 
on  the  computed  undulation  reference  surface,  the  accuracy  of  the  flat- 
earth  approximation  at  these  short  wavelengths  is  readily  apparent. 

In  order  to  illustrate  the  effect  of  survey  track  direction  and 

the  utility  of  these  frequency  domain  operations,  these  procedures  were 

applied  to  observed  free-air  gravity  data  collected  in  a small  test 

area  from  a detailed  shipboard  survey  operation.  Figure  38  is  a 
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Figur*  30.  Ref»r«nc«  gravity  field— 81  X 81  minute  grid  - contour  intervol  {5  mgoU). 
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Figure  31.  Two-dimensional  spectrum  of  model  field  (grovity). 
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Figur*  32.  Actual  two-dim«n$iooal  ETA  spectrum  from  modal. 
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Figure  33.  Actual  two-dimensional  XI  spectrum  from  model. 
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Figuro  34.  Actual  two-dim«nsionol  undulation  ipoctrum  from  modol. 
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Figure  36.  XI  spectrum  computed  with  flat-earth  frequency  domain  operator. 


83 


0.4- 


Figure  37.  Undulation  spectrum  computed  with  flat-earth  frequency  domain  operator. 
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Figure  38.  Gravity  base  data— 21  X 21  point  grid,  1 data  intervol=1  NM, 

contour  interval  (5  mgalt). 
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contour  chart  of  these  gridded  data  at  a grid  interval  of  one  nautical 
mile.  These  data  are  considered  to  represent  a type  area  from  a large 
homogeneous  province  developed  in  the  manner  outlined  in  Figure  22. 

Figure  39  is  a contour  chart  of  the  two-dimensional  amplitude  spectrum 
of  these  data  after  appropriate  prewhitening  and  correction.  This 
figure  indicates  that,  except  for  a small  amount  of  short  wavelength 
energy,  nearly  all  of  the  strong  lineation  shown  in  Figure  38  is  con- 
tained in  wavelengths  much  longer  than  the  dimension  of  the  base  data. 

As  a consequence,  most  of  the  lineation  has  been  removed  by  the  pre- 
whitening process. 

Applying  equation  (6)  to  the  real  and  imaginary  parts  of  this  two- 
dimensional  transform,  with  S^(x,y)  oriented  to  model  north-south  survey 
tracks  at  various  spacings,  results  in  an  estimate  of  the  two-dimensional 
spectral  content  of  the  sampling  error  for  the  measured  gravity  data. 
Figures  AO  and  A1  are  contour  charts  of  the  amplitude  spectrum  of  this 
sampling  error  for  north-south  tracks  spaced  16  nautical  miles  apart 
and  5.3  nautical  miles  apart.  The  lack  of  long  wavelength  lineations 
in  the  error  spectrum  at  the  16  mile  spacing  is  a consequence  of  the 
Isotropic  character  of  the  low-frequency  energy  in  the  base  data.  As 
would  be  expected,  the  error  of  omission  indicated  by  the  contour 
levels  1,  5,  and  10  is  nearly  identical  to  the  energy  in  the  base  data 
at  these  frequencies. 

The  application  of  equations  (1A)-(16)  to  these  estimates  of  the 

gravity  error  spectrum  results  in  straightforward  propagation  of  the 

gravity  error  through  equations  (12)  and  (13)  without  the  need  for 

performing  the  two-dimensional  convolution  operations.  Figures  A2  and 

A3  are  the  resulting  estimates  of  the  sampling  error  spectrum  for 
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Figure  39.  Spectrum  of  gravity  base  data. 
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Fiflor#  40,  Sampling  error  jpectrum  {gravity)  N-S  frocki— 1 6 NM  spacing. 
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Rfluro  41.  Sampling  •rror  ipoctrum  (gravity)  N-S  trocki— 5.3  NM  ipodng. 
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Figure  42.  Sampling  error  spectrum  (XI)  N-S  tracks— 16  NM  spacing. 
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Figur*  43.  Sampling  error  spectrum  (XI)  N-S  tracks— 5.3  NM  spacing. 
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the  along- track  deflection  component.  In  this  case,  for  north-south 
tracks,  this  is  the  XI  component  of  deflection.  Since  the  frequency 
response  of  the  XI  operator  attenuates  all  spectral  components  except 
those  in  the  north-south  direction,  these  error  spectra  show  a pro- 
nounced lineation  at  all  frequencies.  As  was  the  case  for  the  gravity 
error  spectra,  the  error  of  omission  is  nearly  the  same  for  both  track 
spaclngs.  Figures  4A  and  45  are  the  error  spectra  for  the  cross-track 
(ETA)  component.  As  was  the  case  with  the  XI  operator,  the  frequency 
response  of  the  ETA  operator  generates  a strongly  lineated  error 
spectrum  at  both  track  spaclngs.  Figures  46  and  47  show  the  error 
spectra  produced  by  applying  the  frequency  response  of  the  geoid  undula- 
tion operator.  Since  this  operator  exhibits  circular  symmetry,  the 
undulation  error  spectrum  is  simply  an  attenuated  version  of  the  gravity 
error  spectrum  at  each  of  the  sample  track  spaclngs. 

Figures  48-55  are  similar  contours  of  the  estimated  error  spectra 
for  the  case  when  S^(x,y)  is  oriented  to  model  east-west  survey  tracks 
at  various  spaclngs.  For  this  track  orientation,  the  along-track  deflec- 
tion component  is  the  ETA  component,  and  the  cross-track  component  is 
the  XI  component.  In  this  area  of  high-frequency  north-south  llneatlons 
of  the  gravity  field,  the  reduction  of  the  sampling  error  occurring  when 
the  survey  tracks  are  oriented  normal  to  the  trend  is  readily  apparent. 

As  in  the  case  for  north-south  tracks,  the  effect  of  the  frequency 
responses  of  the  various  operators  is  readily  apparent. 

The  result  of  applying  equation  (8)  to  each  of  these  error  spectra 
is  shown  in  Figures  56  and  57.  As  in  the  case  of  the  error  spectra,  the 
influence  of  track  direction  on  the  amplitude  of  the  sampling  error  is  quite 
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Figure  44.  Sampling  eiTor  spectrum  (ETA)  N-S  trock»— 16  NM  tpoong. 
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Figor#  45.  Sampling  •rror  tpoctnim  (ETA)  N-S  track*— 5.3  NM  spacing. 
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Figure  46.  Sampling  error  spectrum  (undulation)  N*S  tracks— 16  NM  spacing. 
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Fiflwr*  47.  Sotnpling  •rror  »p«ctTom  (undulation)  N-S  tracks— 5.3  NM  spacing. 
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Figur*  48.  Sampling  arror  tpactrum  (gravity)  E-W  tracks— 1 6 NM  spacing. 
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Figure  49.  Sampling  •rror  spectrum  (gravity)  E-W  tracks— 5.3  NM  spodng. 
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Figure  50.  Sampling  error  spectrum  (XI)  E-W  tracks — 16NM  spacing. 
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Figure  51 . Sampling  error  spectrum  (XI)  E-W  tracks— 5.3  NM  spacing. 
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Figure  52.  Sampling  error  ipectrum  (ETA)  E-W  tracks— 16  NM  spodng. 
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Figure  53.  Sampling  error  spectrum  (ETA)  E-W  tracks— 5.3  NM  spacing. 
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Figure  55.  Sampling  error  spectrum  (undulotioo)  E-W  tracks — 5.3  NM  spacing. 
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Figura  56.  Sampling  arror  vs.  track  spacing  far  gravity  and  undulation. 
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Figure  57.  Sampling  error  vi.  track  spocing  for  XI  orxl  ETA. 
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pronounced-  In  addition,  it  should  be  noted  that,  in  general,  the 
along-track  component  of  vertical  deflection  contains  less  error  than 
the  cross-track  component.  Reference  to  Figures  11  and  29  indicates 
that  this  numerical  result  is  entirely  consistent  with  that  predicted 
by  the  theory- 


Survey  Design  for  Oceanic  Sound  Speed 

The  procedure  which  has  been  developed  for  estimating  the  bound- 
aries of  relatively  homogeneous  oceanic  sound  speed  provinces  from 
existing  historical  data  is  outlined  in  Figure  23.  Although  oceanic 
sound  speed  is,  in  reality,  a four-dimensional  function,  ie.,  a func- 
tion of  latitude,  longitude,  de-pth,  and  time,  the  distribution  of 
existing  data,  and  the  cost  of  large  scale  shipboard  survey  operations 
precludes  the  simultaneous  consideration  of  all  four  variables  on 
anything  but  a small  test  basis. 

The  basic  approach  to  the  design  of  sound  speed  surveys  consists 
of  applying  the  procedure  for  the  design  of  track-type  surveys  to 
detailed  measurements  made  along  a short  survey  line  within  each  homo- 
geneous province.  In  this  operation,  the  survey  ’’tracks**  are  vertical 
profiles  of  sound  speed  as  a function  of  depth.  In  order  to  test  this 
survey  design  procedure,  a series  of  measurements  were  made  within  a 
particular  homogeneous  winter  season  province  which  was  expected  to 
exhibit  a significant  amount  of  high  frequency  temporal  variability. 

To  obtain  the  time  series  test  data,  17  vertical  profiles  of  sound 
speed  were  collected  at  approximately  the  same  geographic  position  in 
as  short  a time  as  the  mechanics  of  the  operation  v;ould  allow.  The 
depth  range  of  the  measurements  was  from  25  meters  to  1500  meters 

covering  a time  span  of  approximately  2900  minutes. 
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Figure  58  is  a contour  chart  of  the  gridded  residual  data  derived 
from  these  measurements.  This  residual  field  is  constructed  by  remov- 
ing the  average  profile  from  each  observed  profile  in  order  to  show  the 
time  variations  in  depth  ranges  in  which  the  sound  speed  possessed  a 
large  vertical  gradient.  The  actual  time  of  each  profile  or  "station" 
is  shown  in  Figure  58.  Immediately  following  the  time  series  measure- 
ments, a spatial  series  was  collected  along  an  east-west  line.  The.e 
data  consist  of  16  profiles  collected  at  a nominal  spacing  of  15  nauti-  . 
cal  miles  covering  a total  distance  of  210  nautical  miles.  Figure  59 
is  a contour  chart  of  the  gridded  residual  sound  speed  field  derived 
from  this  set  of  measurements.  In  both  pf  these  data  sets,  the  grid 
Interval  in  depth  was  25  meters.  For  the  temporal  variation  data,  the 
horizontal  grid  interval  was  25  minutes,  and  for  the  spatial  variation, 
the  grid  was  1.45  nautical  miles  which  is  roughly  equivalent  to  25 
minutes.  Figures  60  and  61  show  the  two  dimensional  amplitude  spectra 
of  these  data. 

During  the  design  phase  of  this,  experimental  survey  operation,  it 
was  anticipated  that  the  temporal  variations  would  exhibit  a sufficient- 
ly stationary  character  over  the  several  days  of  survey  operation  so 
that  the  temporal  variations  could  be  separated  from  the  spatial  varia- 
tion spectrum  to  produce  a relatively  uncontaninated  space  series.  As 
shown  in  Figures  60  and  61,  such  was  not  the  case.  The  spectral  content 
of  the  time  variations  for  frequencies  less  than  0.2  cycles/data  inter- 
val was  significantly  larger  than  that  obtained  for  the  space  series. 

As  a consequence,  stationarity  could  not  be  assumed  and  the  space 
series  was  analyzed  in  its  original  form.  The  result  of  applying  the 

two-dimensional  survey  design  procedure  to  each  of  these  data  sets  is 
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Figure  60.  Spectrum  of  sound  speed  femporol  variations. 
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Figur*  61 . Sp«ctn»m  of  sound  sp««d  spotiol  voriotions. 
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shown  in  Figure  62.  This  figure  shows  the  estimated  sampling  error 
over  the  depth  range  25-1500  meters  as  a function  of  profile  spacing  in 
both  time  and  space. 

These  data  may  also  be  utilized  to  illustrate  the  facility  with 
which  this  sampling  theory  may  be  adapted  to  many  different  situations. 
In  this  example,  if  the  one-dimensional  Fourier  transform  is  computed 
using  data  at  each  grid  point  occurring  at  a common  depth,  contours  of 
the  amplitude  spectrum  as  a function  of  both  depth  and  normalized 
frequency  in  the  horizontal  dimension  may  be  obtained.  Figures  63  and 
6A  Illustrate  the  resulting  spectra.  A straightforward  application  of 
equations  (3)  and  (4)  to  these  one-dimensional  transforms  results  in  a 
contour  chart  of  estimated  sampling  error  as  a function  of  both  depth 
and  horizontal  spacing  as  sho^m  in  Figures  65  and  66. 

By  comparing  the  average  sampling  error  shown  in  Figure  62  with 
that  shown  in  Figure  65  at  various  depths,  the  importance  of  the  addi- 
tional information  obtained  by  the  one- dimensional  analysis  is  readily 
apparent.  For  example,  assume  that  a temporal  variation  survey  is  to 
be  conducted  in  this  area  to  a specified  accuracy  of  1 (m/sec) Figure 
62  indicates  that  the  required  sample  spacing  to  achieve  this  accuracy 
is  approximately  1050  minutes.  However,  Figure  65  indicates  that  this 
sample  spacing  will  be  expected  to  produce  almost  four  times  the  allow- 
able error  at  depths  less  than  200  meters.  Thus,  the  appropriate  survey 
design  for  this  accuracy  specification  would  consist  of  stations  extend- 
ing to  200  meters  and  spaced  300  minutes  apart  with  every  fourth  station 
extending  to  a depth  of  1500  meters. 
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Figure  62.  Sampling  error  vs.  station  spacing  for  oceanic  sound  speed. 
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Figure  63.  Spectrum  of  sound  speed  temporal  variations  vs.  depth. 
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Figure  64.  Spectrum  of  sound  speed  spatial  variations  vs.  depth. 
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Figure  65.  Samplirtg  error  vs.  depth  For  sound  speed  temporal  voriotions. 


117 


OCPTH  (AAETtRS) 


Figure  66.  Sampling  error  vs.  depth  for  sound  speed  spatial  variations. 
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An  Example  of  the  Effect  of  Sampling  Error 
on  Geophysical  Interpretation 

In  some  cases,  the  magnitude  of  the  sampling  error  generated  by  a 
particular  survey  pattern  may  be  quite  small  and  yet  have  a significant 
effect  upon  operations  performed  on  the  survey  data.  This  example  was 
designed  to  illustrate  the  effect  of  a very  small  sampling  error  upon 
the  interpretation  of  the  gravity  anomaly  associated  with  a two- 
dimensional  fault.  The  interpretation  technique  which  was  applied  to 
equally  spaced  digital  data  computed  along  a profile  normal  to  the 
strike  of  the  fault  model  was  developed  by  Davis  (1971) . This  technique 
is  based  on  the  use  of  characteristic  curves  derived  through  two  digital 
band-pass  filters.  The  amplitude  responses  of  these  filters  were 
selected  to  produce  an  efficient  regional-residual  separation  by  mini- 
mizing the  effect  of  the  overlap  of  anomaly  spectrum  and  regional 
spectrum  by  retaining  only  enough  short  wavelength  anomaly  information 
to  allow  a quantitative  interpretation. 

The  parameters  of  the  fault  model  (Figure  2)  selected  for  this 
example  are  as  follows:  depth  = 5 km,  thickness  = 6 km,  density  con- 

trast =0.5  g/cc,  dip  = 60®  SE,  and  strike  = 20®.  The  gravity  anomaly 
generated  by  this  model  was  computed  for  an  80x80  km  area  at  a grid 
spacing  of  1 km  utilizing  the  formulation  given  in  Davis  (1971) . Figure 
67  is  a contour  chart  of  this  model  anomaly.  The  two-dimensional  ampli- 
tude spectrum  of  this  field  computed  via  the  FFT  and  appropriate 
prewhitening  is  shown  in  Figure  68.  The  strong  lineation  of  the  anomaly 
field  from  this  model  is  reflected  as  a lineation  of  the  two-dimensional 
spectrum  rotated  by  90®.  The  low-amplitude  lineations  parallel  to  the 
frequency  axes  are  a consequence  of  a small  roundoff  error  in  computing 


the  original  gridded  data. 
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Figure  67.  Gravity  field  generated  by  a two-dimensional  fault  model 
> contour  Interval  (5  mgals). 
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Figure  68.  Two*dimensionol  spectrum  of  gravity  fault  model. 
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Tn  order  to  illustrate  the  effect  of  sampling  error  on  the  inter- 
pretation of  this  aiDTialy,  equations  (6)-(8)  were  applieo  to  che  real 
and  ir?J  inry  • .s  of  this  tv o-dimensional  spectrum,  with  S^(x,y) 
n;  • ' 0 ao ’pj.  north-south  survey  tracks  at  various  spacings.  Figu  e 

69  shows  the  resulting  estimates  of  mean  square  sampling  error  as  a 
function  of  north-south  track  spacing.  Note  particularly  the  low  level 
of  sampl'  g error  asr.^  t ^ track  spacing  of  t-‘twe  8 a’ 

6.4  km.  Figure  70  is  a contour  chart  of  the  predicted  two-dimensional 
.•mnlitude  spectrum  of  the  sampl_.»g  error  for  a north-south  track  spacing 
of  S.  ')  km.  In  this  case,  the  cross  track  replication  of  the  transform 
resulted  in  a slight  clockwise  rotation  of  the  main  lineation  of  the 
error  spectrum  with  respect  to -the  direction  of  the  lineation  of  the 
anomalous  field.  This  rotation  indicates  that  the  lineation  of  the 
sampling  error  in  the  two-dimensional  space  domain  will  be  at  a slight 
angle  to  the  strike  of  the  model. 

In  order  to  determine  the  effect  of  this  sampling  error  on  inter- 
pretation, a track-type  survey  v»as  simulated  by  selecting  every  sixth 
colimin  of  data  from  the  model  field  shovm  in  Figure  67  and  utilizing  a 
spline  interpolation  (Davis  and  Kontis,  1970)  to  produce  data  on  a one 
km  grid.  A cubic  spline  interpolation  algorithm  which  was  specifically 
designed  to  extract  a profile  from  a set  of  gridded  data  (Vanwyckhouse, 
1973)  was  then  utilized  with  this  interpolated  grid  to  obtain  a sample 
profile  normal  to  the  strike  of  the  model.  Figure  71  is  a comparison 
of  the  true  one-dimensional  amplitude  spectrum  of  a profile  normal  to 
the  strike  of  the  body  with  the  spectrum  of  the  profile  generated  from 
the  simulated  survey  data.  The  70  percent  response  points  of  the  two 
digital  band-pass  filters  utilized  in  the  interpretation  arc  also  shown 
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N-t  TKACK  tfAONG  (KM) 

Figure  69.  AAeon  square  Mm(9ling  error  vs.  N-S 
track  spacing  for  fauH  model. 
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Figure  70.  Sompling  error  spectrum  for  N-S  trock  spocing  of  5.8  km. 
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NOItMAlIZtO  ntEQOENCY  (VW) 

Figure  71.  Comparison  of  true  spectrum  of  a profile  rtormal  to  strike  of  model  vs. 
spectrum  of  profile  constructed  from  N-S  tracks  spaced  at  6.0  km. 
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In  this  figure.  Note  particularly  the  rather  small  magnitude  of  the 
sampling  error  falling  within  the  pass  band  of  the  filters.  Applying 
the  Interpretation  technique  (Davis,  1971)  to  the  model  profile  which 
was  uncontaminated  by  sampling  error  resulted  in  the  following  estimates 
depth  = A.9  km,  thickness  =6.1  km,  dip  = 61®,  and  density 
contrast  = 0.A5  g/cc.  The  interpreted  quantities  derived  ftom  the 
profile  recovered  from  the  simulated  survey  data  were  as  follows: 
depth  = 6.3  km,  thickness  = 2.3  km,  dip  = 30®,  and  density  contrast  = 
1.10  g/cc.  This  result  clearly  indicates  the  advisability  of  estimating 
sampling  error,  especially  in  cases  where  the  short  wavelength  compon- 
ents of  the  survey  data  are  of  primary  importance. 
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SUMMARY  AND  CONCLUSIONS 


The  requirement  for  designing  geophysical  surveys  in  such  a way 
that  the  true  continuous  field  may  be  derived  to  a predetermined  accu- 
racy from  digital  samples  is  basic  to  any  data  collection  operation. 

The  time  required  in  the  collection  and  analysis  of  reconnaissance 
data,  as  well  as  the  detailed  surveys  of  the  small  type  areas,  may  be 
justified  not  only  from  a purely  scientific  standpoint  but,  with  the 
present  cost  of  shipboard  operations  exceeding  $2000  a day,  from  an 
economic  standpoint  as  well. 

Detailed  tests  of  the  mathematical  procedures  applicable  to  the 
design  of  geophysical  surveys  which  have  been  presented  indicate  the 
accuracy  and  efficiency  with  which  the  appropriate  estimates  of 
sampling  error  may  be  derived  by  operating  in  the  two-dimensional 
frequency  domain.  The  importance  of  some  form  of  prewhitening  or 
window  modification  for  computing  accurate  spectral  estimates  from  a 
finite  length  of  data  obtained  from  a "red  noise”  process  has  been 
clearly  demonstrated.  The  raster  sampling  model  consisting  of  a set  of 
parallel  delta  function  ridges  has  been  shown  to  be  an  appropriate 
mathematical  model  for  defining  two-dimensional  track-type  surveys. 

Several  practical  applications  have  been  presented  to  demonstrate 
the  adaptability  of  the  theory  to  a wide  variety  of  survey  operations 
and  the  relative  ease  with  which  estimates  of  mean  square  error  and 
the  two-dimensional  spectral  content  of  this  sampling  error  may  be 
computed.  The  use  of  a so-called  power  fold  approximation,  combined 
with  a priori  knowledge  of  the  orientation  of  major  lineations,  has 
been  developed  for  survey  design  in  near  real-time  using  small  scale 


127 


computers.  The  practical  procedures  which  are  presented  for  designing 
surveys  to  produce  products  of  a specified  accuracy  in  which  the  accu- 
racy requirement  is  placed  not  on  the  survey  data  but  upon  quantities 
computed  f -om  the  data  by  two-dimensional  cross  correlation  or  convolu- 
tion operations  is  of  particular  importance  from  the  standpoint  of 
systems  analysis.  Finally,  an  example  of  the  effect  of  sampling  error 
on  a problem  of  geophysical  interpretation  clearly  illustrates  the  fact 
that,  although  the  magnitude  of  the  sampling  error  may  be  quite  small, 
it  may  have  a significant  effect  upon  operations  performed  on  the 
survey  data. 

No  attempt  has  been  made  to  specify  an  allowable  sampling  error 
for  any  of  the  survey  designs  which  have  been  presented  here.  This 
decision  can  only  be  made  on  a case  by  case  basis  and  will  depend  upon 
the  ultimate  use  of  the  survey  data.  The  test  results  which  have  been 
presented  indicate  clearly  that  this  theory  will,  for  the  first  time, 
provide  the  individual  investigator  with  an  accurate  and  efficient 
technique  for  designing  track-type  geophysical  surveys  to  meet  a pre- 
determined accuracy  specification. 

Future  work  in  this  area  will  be  directed  toward  modifying  the 
procedure  for  designing  surveys  for  vertical  deflections  in  order  to 
remove  the  requirement  for  detailed  type  area  surveys.  This  modifica- 
tion should  provide  a basis  for  real-time  survey  design  through  the  use 
of  reconnaissance  tracks  normal  to  the  survey  pattern  combined  with  an 
assumption  of  short  wavelength  isotropy  to  produce  an  estimate  of  the 
two-dimensional  spectrum  which  will  be  accurate  in  the  cross-track 
direction.  , 
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APPENDIX  A 


THE  DESIGN  OF  PREWUTENING  FILTERS 

An  excellent  technique  for  the  practical  computation  of  the  digital 
weight  function  required  for  one-diraensional  prewhitening  filters  has 
been  developed  by  Martin  (1957)  and  evaluated  by  Linnette  (1961).  For 
general  low-pass  filters,  the  2N+1  filter  weights,  symmetric  about  K ■ 0 
to  Insure  no  phase  shift,  are  computed  from 


Cos  (2TTKH) 

Sin  2ttK  (V  +H)“ 
c 

[l-16H^K^ 

TTK 

K = -N....0....N, 


(A-1) 


where  K is  the  weight  number,  H -is  a parameter  controlling  the  shape  of 
the  amplitude  response,  and  is  the  normalized  cutoff  frequency  in 
cycles/data  interval.  The  parameter  H is  usually  assigned  a value 
between  0.01  and  0.20.  In  general,  for  a specified  value  of  N,  the 
slope  of  the  frequency  response  of  the  filter  and  the  magnitude  of  the 
first  side  lobe  decrease  as  H is  increased.  The  cutoff  frequency  is 
defined  as  that  point  where  the  amplitude  response  leaves  100  percent  in 
the  case  of  a low-pass  filter  or  reaches  zero  in  the  case  of  a high-pass 
filter.  In  order  for  the  filter  to  have  unity  gain,  a correction  must 
be  applied  to  the  weights  generated  by  equation  (A-1)  to  insure  that 
their  sum  is  unity.  This  correction  is 


A 


(A-2) 


At  values  of  K = 0 and  K = ± H/A,  equation  (A-1)  is  of  indeter- 
minate form.  Application  of  L’Hospital’s  rule  to  evaluate  the  limit  as 
K approaches  these  values  yields 
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and 


L “ 2(V  +H)  for  K = 0 
o c 


Sin(2TTK(V  +H)] 

h.  “ for  K = ± H/4. 


The  final  corrected  filter  weights  for  the  general  low-pass  digital 

A 


filter  arv  given  by  W, 


K 


2N+1* 


The  weights  for  the  high-pass  prewhitening  filter  (P  ) are  simply 

constru»'ted  by  subtracting  W from  an  all-pass  filter  defined  as 

K 


K 


[1,  K=0 

10,  K=±l ±N  . 


Thus,  P 


K 


1-Wj^.  K=0 

-W  , K=±l ±N 

K. 


Experience  has  shown  that,  for  one-dimensional  prewhitening  of  a 
gravity  or  magnetic  profile  containing  T digital  values,  the  filter  con- 
trol parameters  given  by  V = — - , H = 0.2,  N = 3 produced  a reason- 

ably  white  spectrum. 

A computational  procedure  for  generating  the  two-dimensional  digi- 
tal weight  function  for  prewhitening  gridded  data  has  been  developed  by 
Lavin  and  Devane  (1970).  The  digital  filter  weights  possess  circular 
symmetry  to  insure  circular  symmetry  in  the  two-dimensional  spatial 
frequency  domain,  and  are  given  by 

aJ  (2Trar)  J (irrAk) 


W(r)  = 


j^_^2]irAKj2  * 


(A-3) 


1/2 

where  r * (x^+y^)  , a = (K  +K  )/2  with  K defined  as  the  normalized 

v«  X O 

cutoff  frequency  and  the  normalized  termination  frequency  in  cycles 
per  data  interval,  a = 4.80965....,  and  AK  * 


130 


Equation  (A-3)  is  undefined  at  r “ 0 and  r = 


g 

2irAK 


. Application  of 


L'Hospital*s  rule  at  these  points  yields  W(o)  = Tia^  and 


" ~~2^  As  in  the  one-dimensional  case,  a correc- 

tion equivalent  to  equation  (A-2)  is  applied  to  insure  unity  gain,  and 
the  high-pass  prewhitening  weights  are  derived  by  subtracting  equation 
(A-3)  from  the  two-dimensional  all-pass  filter.  In  general,  an  accept- 
able prewhitening  filter  for  gridded  "red  noise"  type  processes,  such  as 

gravity  or  magnetic  data,  may  be  generated  by  setting 
1 


K 


r,  AK  = 0.3,  and  terminating  equation  (A-3)  by  setting  the 


c T+O.IT’ 

maximum  value  of  x and  y equal  to  four  data  intervals,  ie. , 


x,y  = -4. . . .0. . . .4. 
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APPENDIX  B 


THE  SHANNON  SAMPLING  THEOREM 


Shannon  (1949)  states  that,  given  equally  spaced  digital 
values  f(mT)  of  a continuous  function  f(t)  which  is  band-limited  to  the 

fir 

normalized  frequency  range  -TT_<aKiT,  and  such  that  I |F(a))ldo)  exists, 

J-Tj 

f(t)  may  be  represented  by 


f(t)  = E f(mT) 

QS-CO 


Sin  Tr(t-mT) 
TT(t-mT) 


for  all  t. 


The  proof  of  this  theorem  consists  of  showing  that 

lira  [f(t)  - E f(mT)  1 = 0 for  all  t. 

-N  Tt(t-mT) 

To  begin  with,  the  inverse  Fourier  transform  of  the  so-called  boxcar 
function,  defined  as 

. fl,  IwUir 

\o,  elsewhere’ 

Is  required.  By  definition, 

b(t)  “ r B(a))e^‘^*'d(o  = ^ 

J <-co 

Thus,  utilizing  the  shifting  theorem  for  Fourier  transforms  (Papoulis, 
1962) , we  have 

Sin  TT(t-mT)  _ 1 
iT(t-mT)  ~ 2ir 


ita(t-mT)  , 
e du). 


cos  CJtdto 


Simrt 

irt 


Then 
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I ^ 

0<  f(t)  -I  f(mT) 

I _N 


Sin  Tr(t-tnT) 
7r(t-mT) 


i_r 

J ^ 

' -TT 


F(u)e‘“'<to  - It  S 

2"  -N 


r 

“’-ir 


ito(t-mT)  , 
e ' ''do) 


2tt 


i:. 


N 


F(w)  - Z f(mX)e  " |dw 
-N 


-iOjraTl 

J' 


(B-l) 


Since  F(o))  is  defined  in  the  interval  it  may  be  expanded  in  a 


^ I'p 

Fourier  series  F(w)=  E c e where  c 

m m 

.00 

since  f(t)  is  bandlimited , 

f(t) 


. io)m, 

F(o))e  do).  But, 


2tt  I 

''-7T 


F(w)  e^^^do). 


60 


c “ f(mT) 
m 


.Lf 

J-F 


. iwmT, 

F (co)  e dw , 


and 


F(o))=  E f(mT)e 

IQS— OO 


-iwmT 


(B-2) 


Applying  the  Schwarz  inequality  (Kaplan,  1952),  given  by 

1/2 


r 

r f,  ,2  i 

1/2 

J f(x)g(x)dx 

< 

J If (x) 1 dx 

J|g(x)|  dx^ 

to  equation  (B-l)  yields 


1 iwtri.  . ? _ -iurafl^  J 

■^r  e F((o)  - E f(nT)e  du 

L m=-N  -J  ' 

-■kff  ^ f(iuT)e"^““^NJ 

'•■'-IT  i '••'-TT  m=-N  ^ 


-iumTi*  11/2  . _ 

e I do)|  = A'B, 
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where 


1 ff^  I Iwti* 


1/2 


and 


(f-n  N 

|f(o))  - E f(mT) 

U_Ti  m=-N 


-IWraTi ^ . 
e do) 


1/2 


Since  |e^^^|  = 1,  then  A = l/ZlH,  and  utilizing  equation  (2),  it  is 

clear  that,  as  N->*,  thus 


£(t)  = lim  E f(mT)  for  all  t. 

IFKo  m=-N 
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